proper linking for libHLTrec.so
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerNxN.cxx
CommitLineData
ee08edde 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// This class derives from AliEMCALClustrerizer but also keeps the API of AliEMCALClusterizerv1
17// Algorithm:
18// 1. peek the most energetic cell
19// 2. assign it as a center of the cluster and add cells surrounding it: 3x3, 5x5...
20// 3. remove the cells contributing to the cluster
21// 4. start from 1 for the remaining clusters
22// 5. cluster splitting (not implemented yet) - use the shape analysis to resolve the energy sharing
23// - for high energy clusters check the surrounding of the 3x3 clusters for extra energy
24// (merge 3x3 clusters and resolve the internal energy sharing - case for 2 clusters merged)
25// Use Case:
26// root [0] AliEMCALClusterizerNxN * cl = new AliEMCALClusterizerNxN("galice.root")
27// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
28// //reads gAlice from header file "..."
29// root [1] cl->ExecuteTask()
30// //finds RecPoints in all events stored in galice.root
31// root [2] cl->SetDigitsBranch("digits2")
32// //sets another title for Digitis (input) branch
33// root [3] cl->SetRecPointsBranch("recp2")
34// //sets another title four output branches
35// root [4] cl->SetTowerLocalMaxCut(0.03)
36// //set clusterization parameters
37// root [5] cl->ExecuteTask("deb all time")
38// //once more finds RecPoints options are
39// // deb - print number of found rec points
40// // deb all - print number of found RecPoints and some their characteristics
41// // time - print benchmarking results
42
43// --- ROOT system ---
44#include <TMath.h>
45#include <TMinuit.h>
46#include <TTree.h>
47#include <TBenchmark.h>
48#include <TBrowser.h>
49#include <TROOT.h>
50#include <TClonesArray.h>
51
52// --- Standard library ---
53#include <cassert>
54
55// --- AliRoot header files ---
56#include "AliLog.h"
57#include "AliEMCALClusterizerNxN.h"
58#include "AliEMCALRecPoint.h"
59#include "AliEMCALDigit.h"
60#include "AliEMCALGeometry.h"
61#include "AliCaloCalibPedestal.h"
62#include "AliEMCALCalibData.h"
63#include "AliESDCaloCluster.h"
65bec413 64#include "AliEMCALUnfolding.h"
ee08edde 65
66ClassImp(AliEMCALClusterizerNxN)
67
ee08edde 68//____________________________________________________________________________
69AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
303ff65c 70 : AliEMCALClusterizer(), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
ee08edde 71{
72 // ctor with the indication of the file where header Tree and digits Tree are stored
73}
74
75//____________________________________________________________________________
76AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
303ff65c 77 : AliEMCALClusterizer(geometry), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
ee08edde 78{
79 // ctor with the indication of the file where header Tree and digits Tree are stored
80 // use this contructor to avoid usage of Init() which uses runloader
81 // change needed by HLT - MP
ee08edde 82}
83
84//____________________________________________________________________________
85AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
303ff65c 86: AliEMCALClusterizer(geometry, calib, caloped), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
ee08edde 87{
8a84b948 88 // ctor, geometry and calibration are initialized elsewhere.
ee08edde 89}
90
ee08edde 91//____________________________________________________________________________
8a84b948 92AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
ee08edde 93{
94 // dtor
95}
96
ee08edde 97//____________________________________________________________________________
98void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
99{
100 // Steering method to perform clusterization for the current event
101 // in AliEMCALLoader
51e4198e 102
bd0552e8 103 if (strstr(option,"tim"))
ee08edde 104 gBenchmark->Start("EMCALClusterizer");
105
bd0552e8 106 if (strstr(option,"print"))
107 Print("");
51e4198e 108
ee08edde 109 //Get calibration parameters from file or digitizer default values.
bd0552e8 110 GetCalibrationParameters();
51e4198e 111
ee08edde 112 //Get dead channel map from file or digitizer default values.
bd0552e8 113 GetCaloCalibPedestal();
ee08edde 114
bd0552e8 115 MakeClusters(); //only the real clusters
51e4198e 116
bd0552e8 117 if (fToUnfold) {
65bec413 118 fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
119 fClusterUnfolding->MakeUnfolding();
120 }
51e4198e 121
65bec413 122 //Evaluate position, dispersion and other RecPoint properties for EC section
bd0552e8 123 for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
51e4198e 124 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
bd0552e8 125 if (rp) {
294553d1 126 rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
51e4198e 127 AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
ee08edde 128 //For each rec.point set the distance to the nearest bad crystal
bd0552e8 129 if (fCaloPed)
130 rp->EvalDistanceToBadChannels(fCaloPed);
ee08edde 131 }
51e4198e 132 }
ee08edde 133
bd0552e8 134 fRecPoints->Sort();
51e4198e 135
bd0552e8 136 for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
137 AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
138 if (rp) {
139 rp->SetIndexInList(index);
7e1d9a9b 140 }
141 else AliFatal("RecPoint NULL!!");
51e4198e 142 }
ee08edde 143
bd0552e8 144 if (fTreeR)
145 fTreeR->Fill();
ee08edde 146
bd0552e8 147 if (strstr(option,"deb") || strstr(option,"all"))
148 PrintRecPoints(option);
51e4198e 149
ee08edde 150 AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
51e4198e 151
bd0552e8 152 if (strstr(option,"tim")) {
ee08edde 153 gBenchmark->Stop("EMCALClusterizer");
154 printf("Exec took %f seconds for Clusterizing",
51e4198e 155 gBenchmark->GetCpuTime("EMCALClusterizer"));
ee08edde 156 }
157}
158
159//____________________________________________________________________________
160Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
161{
c526514b 162 // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
163 // = 1 are neighbour
164 // = 2 is in different SM; continue searching
165 // In case it is in different SM, but same phi rack, check if neigbours at eta=0
166 // neighbours are defined as digits having at least a common side
167 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
168 // which is compared to a digit (d2) not yet in a cluster
169
303ff65c 170 if (fEnergyGrad) { //false by default
171 if (d2->GetCalibAmp()>d1->GetCalibAmp())
172 return 3; // energy of neighboring cell should be smaller in order to become a neighbor
173 }
174
bd0552e8 175 Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
176 Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
177 Int_t rowdiff=0, coldiff=0;
c526514b 178
179 shared = kFALSE;
180
181 fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
182 fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
183 fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
184 fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
185
186 //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
bd0552e8 187 if (nSupMod1 != nSupMod2 )
c526514b 188 {
189 //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
190 Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
191 Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
192
bd0552e8 193 if (!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
c526514b 194
195 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
196 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
bd0552e8 197 if (nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
c526514b 198 else ieta2+=AliEMCALGeoParams::fgkEMCALCols;
199
200 shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
201
202 }//Different SM, same phi
203
204 rowdiff = TMath::Abs(iphi1 - iphi2);
bd0552e8 205 coldiff = TMath::Abs(ieta1 - ieta2);
c526514b 206
207 // neighbours +-1 in col and row
8a84b948 208 if ( TMath::Abs(coldiff) <= fNColDiff && TMath::Abs(rowdiff) <= fNRowDiff)
c526514b 209 {
210
211 AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
212 d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
213
214 return 1;
215 }//Neighbours
216 else
217 {
218 AliDebug(9, Form("NOT AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
219 d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
220 shared = kFALSE;
bd0552e8 221 return 2;
c526514b 222 }//Not neighbours
ee08edde 223}
224
225//____________________________________________________________________________
226void AliEMCALClusterizerNxN::MakeClusters()
227{
8a84b948 228 // Make clusters
51e4198e 229
8a84b948 230 if (fGeom==0)
231 AliFatal("Did not get geometry from EMCALLoader");
51e4198e 232
bd0552e8 233 fNumberOfECAClusters = 0;
234 fRecPoints->Delete();
51e4198e 235
783153ff 236 // Set up TObjArray with pointers to digits to work on, calibrate digits
ee08edde 237 TObjArray digitsC;
238 TIter nextdigit(fDigitsArr);
239 AliEMCALDigit *digit = 0;
bd0552e8 240 while ( (digit = static_cast<AliEMCALDigit*>(nextdigit())) ) {
783153ff 241 Float_t dEnergyCalibrated = digit->GetAmplitude();
242 Float_t time = digit->GetTime();
243 Calibrate(dEnergyCalibrated,time ,digit->GetId());
8a84b948 244 digit->SetCalibAmp(dEnergyCalibrated);
783153ff 245 digit->SetTime(time);
ee08edde 246 digitsC.AddLast(digit);
247 }
51e4198e 248
ee08edde 249 TIter nextdigitC(&digitsC);
51e4198e 250
ee08edde 251 AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f)\n",
51e4198e 252 fDigitsArr->GetEntries(),fMinECut));
253
ee08edde 254 Bool_t bDone = kFALSE;
255 while ( bDone != kTRUE )
51e4198e 256 {
257 //first sort the digits:
258 Int_t iMaxEnergyDigit = -1;
259 Float_t dMaxEnergyDigit = -1;
260 AliEMCALDigit *pMaxEnergyDigit = 0;
261 nextdigitC.Reset();
8a84b948 262 while ( (digit = static_cast<AliEMCALDigit *>(nextdigitC())) )
51e4198e 263 { // scan over the list of digitsC
8a84b948 264 Float_t dEnergyCalibrated = digit->GetCalibAmp();
959d5ca4 265 Float_t time = digit->GetTime();
266 if (fGeom->CheckAbsCellId(digit->GetId()) &&
267 dEnergyCalibrated > fMinECut &&
268 time < fTimeMax &&
269 time > fTimeMin ) // no threshold by default!
270 { // needs to be set in OCDB!
8a84b948 271 if (dEnergyCalibrated > dMaxEnergyDigit)
51e4198e 272 {
273 dMaxEnergyDigit = dEnergyCalibrated;
274 iMaxEnergyDigit = digit->GetId();
275 pMaxEnergyDigit = digit;
276 }
8a84b948 277 }
51e4198e 278 }
51e4198e 279 if (iMaxEnergyDigit < 0 || digitsC.GetEntries() <= 0)
280 {
281 bDone = kTRUE;
282 continue;
283 }
284
8a84b948 285 AliDebug (2, Form("Max digit found: %1.5f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
51e4198e 286
287 // keep the candidate digits in a list
288 TList clusterDigitList;
289 clusterDigitList.SetOwner(kFALSE);
290 clusterDigitList.AddLast(pMaxEnergyDigit);
291
292 Double_t clusterCandidateEnergy = dMaxEnergyDigit;
293
8a84b948 294 // now loop over the rest of the digits and cluster into NxN cluster
51e4198e 295 // we do not actually cluster yet: we keep them in the list clusterDigitList
296 nextdigitC.Reset();
bd0552e8 297 while ( (digit = static_cast<AliEMCALDigit *>(nextdigitC())) )
51e4198e 298 { // scan over the list of digitsC
299 if (digit == pMaxEnergyDigit) continue;
8a84b948 300 Float_t dEnergyCalibrated = digit->GetCalibAmp();
51e4198e 301 AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
bd0552e8 302 if (fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0 )
8a84b948 303 {
304 Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR?
bd0552e8 305 if (TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR?
8a84b948 306 Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
307 if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!!
51e4198e 308 {
bd0552e8 309 clusterDigitList.AddLast(digit);
51e4198e 310 clusterCandidateEnergy += dEnergyCalibrated;
311 }
8a84b948 312 }
51e4198e 313 }// loop over the next digits
314
315 // start a cluster here only if a cluster energy is larger than clustering threshold
51e4198e 316 AliDebug(5, Form("Clusterization threshold is %f MeV", fECAClusteringThreshold));
317 if (clusterCandidateEnergy > fECAClusteringThreshold)
318 {
bd0552e8 319 if (fNumberOfECAClusters >= fRecPoints->GetSize())
320 fRecPoints->Expand(2*fNumberOfECAClusters+1);
ee08edde 321
c828bc97 322 (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
323 AliEMCALRecPoint *recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
324 // AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
325 // fRecPoints->AddAt(recPoint, fNumberOfECAClusters);
326 // recPoint = static_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters));
bd0552e8 327 if (recPoint) {
328 fNumberOfECAClusters++;
51e4198e 329 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
330
331 AliDebug(9, Form("Number of cells per cluster (max is 9!): %d", clusterDigitList.GetEntries()));
332 for (Int_t idig = 0; idig < clusterDigitList.GetEntries(); idig++)
333 {
51e4198e 334 digit = (AliEMCALDigit*)clusterDigitList.At(idig);
51e4198e 335 AliDebug(5, Form(" Adding digit %d", digit->GetId()));
336 // note: this way the sharing info is lost!
783153ff 337 recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
51e4198e 338 digitsC.Remove(digit);
339 }
340 }// recpoint
341 }
342 else
343 {
344 // we do not want to start clustering in the same spot!
345 // but in this case we may NOT reuse this seed for another cluster!
346 // need a better bookeeping?
347 digitsC.Remove(pMaxEnergyDigit);
348 }
349
350 AliDebug (2, Form("Number of digits left: %d", digitsC.GetEntries()));
351 } // while ! done
352
ee08edde 353 AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
354}