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