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 : //-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
17 : //-- Gustavo Conesa (LPSC-Grenoble), move common clusterizer functionalities to mother class
18 : //////////////////////////////////////////////////////////////////////////////
19 : // Clusterization class. Performs clusterization (collects neighbouring active cells) and
20 : // unfolds the clusters having several local maxima.
21 : // Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
22 : // EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all
23 : // parameters including input digits branch title, thresholds etc.)
24 : //
25 :
26 : // --- ROOT system ---
27 :
28 : #include <TFile.h>
29 : #include <TMath.h>
30 : #include <TMinuit.h>
31 : #include <TTree.h>
32 : #include <TBenchmark.h>
33 : #include <TBrowser.h>
34 : #include <TROOT.h>
35 : #include <TList.h>
36 : #include <TClonesArray.h>
37 :
38 : // --- Standard library ---
39 : #include <cassert>
40 :
41 : // --- AliRoot header files ---
42 : #include "AliLog.h"
43 : #include "AliEMCALClusterizerv1.h"
44 : #include "AliEMCALRecPoint.h"
45 : #include "AliEMCALDigit.h"
46 : #include "AliEMCALGeometry.h"
47 : #include "AliCaloCalibPedestal.h"
48 : #include "AliEMCALCalibData.h"
49 : #include "AliEMCALCalibTime.h"
50 : #include "AliESDCaloCluster.h"
51 : #include "AliEMCALUnfolding.h"
52 :
53 42 : ClassImp(AliEMCALClusterizerv1)
54 :
55 : //____________________________________________________________________________
56 0 : AliEMCALClusterizerv1::AliEMCALClusterizerv1(): AliEMCALClusterizer()
57 0 : {
58 : // ctor with the indication of the file where header Tree and digits Tree are stored
59 :
60 0 : Init();
61 0 : }
62 :
63 : //____________________________________________________________________________
64 : AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry)
65 0 : : AliEMCALClusterizer(geometry)
66 0 : {
67 : // ctor with the indication of the file where header Tree and digits Tree are stored
68 : // use this contructor to avoid usage of Init() which uses runloader
69 : // change needed by HLT - MP
70 0 : }
71 :
72 : //____________________________________________________________________________
73 : AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry,
74 : AliEMCALCalibData * calib,
75 : AliEMCALCalibTime * calibt,
76 : AliCaloCalibPedestal * caloped)
77 2 : : AliEMCALClusterizer(geometry, calib, calibt, caloped)
78 10 : {
79 : // ctor, geometry and calibration are initialized elsewhere.
80 4 : }
81 :
82 : //____________________________________________________________________________
83 : AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
84 8 : {
85 : // dtor
86 8 : }
87 :
88 : //____________________________________________________________________________
89 : void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
90 : {
91 : // Steering method to perform clusterization for the current event
92 : // in AliEMCALLoader
93 :
94 16 : if(strstr(option,"tim"))
95 0 : gBenchmark->Start("EMCALClusterizer");
96 :
97 8 : if(strstr(option,"print"))
98 0 : Print("");
99 :
100 : //Get calibration parameters from file or digitizer default values.
101 8 : GetCalibrationParameters();
102 :
103 : //Get dead channel map from file or digitizer default values.
104 8 : GetCaloCalibPedestal();
105 :
106 8 : fNumberOfECAClusters = 0;
107 :
108 8 : MakeClusters(); //only the real clusters
109 :
110 8 : if(fToUnfold)
111 : {
112 0 : fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
113 0 : fClusterUnfolding->MakeUnfolding();
114 0 : }
115 :
116 : //Evaluate position, dispersion and other RecPoint properties for EC section
117 : Int_t index;
118 82 : for(index = 0; index < fRecPoints->GetEntries(); index++)
119 : {
120 99 : AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
121 :
122 33 : if(rp)
123 : {
124 33 : rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
125 :
126 : // Calculate the number of local maxima in cluster
127 : // Do not do it for unfolded clusters
128 33 : if(!fToUnfold)
129 : {
130 33 : Int_t nMax = rp->GetNumberOfLocalMax(rp->GetMultiplicity(),fECALocMaxCut,fDigitsArr) ;
131 :
132 33 : rp->SetNExMax(nMax);
133 33 : }
134 :
135 : // For each rec.point set the distance to the nearest bad crystal
136 33 : if (fCaloPed)
137 33 : rp->EvalDistanceToBadChannels(fCaloPed);
138 : }
139 0 : else AliFatal("Null rec point in list!");
140 : }
141 :
142 8 : fRecPoints->Sort();
143 :
144 82 : for(index = 0; index < fRecPoints->GetEntries(); index++) {
145 99 : AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
146 33 : if(rp){
147 33 : rp->SetIndexInList(index);
148 33 : rp->Print();
149 33 : }
150 0 : else AliFatal("Null rec point in list!");
151 : }
152 :
153 8 : if (fTreeR)
154 8 : fTreeR->Fill();
155 :
156 16 : if(strstr(option,"deb") || strstr(option,"all"))
157 0 : PrintRecPoints(option);
158 :
159 24 : AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
160 :
161 8 : if(strstr(option,"tim")){
162 0 : gBenchmark->Stop("EMCALClusterizer");
163 0 : printf("Exec took %f seconds for Clusterizing",
164 0 : gBenchmark->GetCpuTime("EMCALClusterizer"));
165 0 : }
166 8 : }
167 :
168 : //____________________________________________________________________________
169 : Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
170 : {
171 : // Gives the neighbourness of two digits = 0 are not neighbour; continue searching
172 : // = 1 are neighbour
173 : // = 2 is in different SM; continue searching
174 : // In case it is in different SM, but same phi rack, check if neigbours at eta=0
175 : // neighbours are defined as digits having at least a common side
176 : // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
177 : // which is compared to a digit (d2) not yet in a cluster
178 :
179 1610 : Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
180 805 : Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
181 :
182 805 : shared = kFALSE;
183 :
184 805 : fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
185 805 : fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
186 805 : fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
187 805 : fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
188 :
189 : //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
190 805 : if (nSupMod1 != nSupMod2 ) {
191 : //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
192 573 : Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
193 573 : Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
194 :
195 1068 : if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
196 :
197 : // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
198 : // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
199 115 : if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
200 41 : else ieta2+=AliEMCALGeoParams::fgkEMCALCols;
201 :
202 78 : shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
203 78 : } //Different SM, same phi
204 :
205 310 : Int_t rowdiff = TMath::Abs(iphi1 - iphi2);
206 310 : Int_t coldiff = TMath::Abs(ieta1 - ieta2);
207 :
208 : // neighbours with at least common side; May 11, 2007
209 710 : if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) {
210 : //Diagonal?
211 : //if ((coldiff==0 && TMath::Abs(rowdiff==1)) || (rowdiff==0 && TMath::Abs(coldiff==1)) || (TMath::Abs(rowdiff)==1 && TMath::Abs(coldiff==1))) rv = 1;
212 :
213 64 : if (gDebug == 2)
214 0 : printf("AliEMCALClusterizerv1::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
215 0 : d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared);
216 64 : return 1;
217 : } //Neighbours
218 : else {
219 246 : shared = kFALSE;
220 246 : return 2;
221 : } //Not neighbours
222 805 : }
223 :
224 : //____________________________________________________________________________
225 : void AliEMCALClusterizerv1::MakeClusters()
226 : {
227 : // Steering method to construct the clusters stored in a list of Reconstructed Points
228 : // A cluster is defined as a list of neighbour digits
229 : // Mar 03, 2007 by PAI
230 :
231 16 : if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
232 :
233 8 : fRecPoints->Delete();
234 :
235 : // Set up TObjArray with pointers to digits to work on calibrated digits
236 8 : TObjArray *digitsC = new TObjArray();
237 : AliEMCALDigit *digit;
238 8 : Float_t dEnergyCalibrated = 0.0, ehs = 0.0, time = 0.0;
239 8 : TIter nextdigit(fDigitsArr);
240 492 : while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // calibrate and clean up digits
241 115 : dEnergyCalibrated = digit->GetAmplitude();
242 115 : time = digit->GetTime();
243 115 : Calibrate(dEnergyCalibrated,time,digit->GetId());
244 115 : digit->SetCalibAmp(dEnergyCalibrated);
245 115 : digit->SetTime(time);
246 :
247 339 : if ( dEnergyCalibrated < fMinECut || time > fTimeMax || time < fTimeMin ){
248 : continue;
249 : }
250 224 : else if (!fGeom->CheckAbsCellId(digit->GetId()))
251 : continue;
252 : else{
253 112 : ehs += dEnergyCalibrated;
254 112 : digitsC->AddLast(digit);
255 : }
256 : }
257 :
258 40 : AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %f\n",
259 : fDigitsArr->GetEntries(),fMinECut,ehs));
260 :
261 8 : TIter nextdigitC(digitsC);
262 329 : while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
263 122 : TArrayI clusterECAdigitslist(fDigitsArr->GetEntries());
264 61 : dEnergyCalibrated = digit->GetCalibAmp();
265 61 : time = digit->GetTime();
266 183 : if(fGeom->CheckAbsCellId(digit->GetId()) && ( dEnergyCalibrated > fECAClusteringThreshold ) ){
267 : // start a new Tower RecPoint
268 66 : if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1);
269 :
270 66 : AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
271 33 : fRecPoints->AddAt(recPoint, fNumberOfECAClusters);
272 132 : recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters));
273 33 : if (recPoint) {
274 33 : fNumberOfECAClusters++;
275 :
276 33 : recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
277 33 : recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
278 33 : TObjArray clusterDigits;
279 33 : clusterDigits.AddLast(digit);
280 33 : digitsC->Remove(digit);
281 :
282 165 : AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(), dEnergyCalibrated, fECAClusteringThreshold)); //Time or TimeR?
283 :
284 : // Grow cluster by finding neighbours
285 33 : TIter nextClusterDigit(&clusterDigits);
286 :
287 584 : while ( (digit = dynamic_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster
288 97 : TIter nextdigitN(digitsC);
289 : AliEMCALDigit *digitN = 0; // digi neighbor
290 1901 : while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours
291 : //Do not add digits with too different time
292 805 : Bool_t shared = kFALSE;//cluster shared by 2 SuperModules?
293 805 : if(TMath::Abs(time - digitN->GetTime()) > fTimeCut ) continue; //Time or TimeR?
294 1610 : if (AreNeighbours(digit, digitN, shared)==1) { // call (digit,digitN) in THAT order !!!!!
295 64 : recPoint->AddDigit(*digitN, digitN->GetCalibAmp(), shared);
296 64 : clusterDigits.AddLast(digitN);
297 64 : digitsC->Remove(digitN);
298 : } // if(ineb==1)
299 1610 : } // scan over digits
300 97 : } // scan over digits already in cluster
301 :
302 165 : AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy()));
303 33 : }//recpoint
304 0 : else AliFatal("Null recpoint in array!");
305 33 : } // If seed found
306 61 : } // while digit
307 :
308 16 : delete digitsC;
309 :
310 40 : AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
311 8 : }
|
