1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // This class derives from AliEMCALClustrerizer but also keeps the API of AliEMCALClusterizerv1
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)
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
43 // --- ROOT system ---
47 #include <TBenchmark.h>
50 #include <TClonesArray.h>
52 // --- Standard library ---
55 // --- AliRoot header files ---
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"
65 ClassImp(AliEMCALClusterizerNxN)
67 Bool_t AliEMCALClusterizerNxN::fgkIsInputCalibrated = kFALSE;
69 //____________________________________________________________________________
70 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
71 : AliEMCALClusterizer()
73 // ctor with the indication of the file where header Tree and digits Tree are stored
76 //____________________________________________________________________________
77 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
78 : AliEMCALClusterizer(geometry)
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
86 //____________________________________________________________________________
87 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
88 : AliEMCALClusterizer(geometry, calib, caloped)
90 // ctor, geometry and calibration are initialized elsewhere.
95 //____________________________________________________________________________
96 AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
102 //____________________________________________________________________________
103 void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
105 // Steering method to perform clusterization for the current event
108 if(strstr(option,"tim"))
109 gBenchmark->Start("EMCALClusterizer");
111 if(strstr(option,"print"))
114 //Get calibration parameters from file or digitizer default values.
115 GetCalibrationParameters() ;
117 //Get dead channel map from file or digitizer default values.
118 GetCaloCalibPedestal() ;
120 fNumberOfECAClusters = 0;
122 MakeClusters() ; //only the real clusters
129 //Evaluate position, dispersion and other RecPoint properties for EC section
130 for(index = 0; index < fRecPoints->GetEntries(); index++)
132 dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalAll(fECAW0,fDigitsArr) ;
133 AliDebug(5, Form("MAX INDEX %d ", dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->GetMaximalEnergyIndex()));
134 //For each rec.point set the distance to the nearest bad crystal
135 dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalDistanceToBadChannels(fCaloPed);
140 for(index = 0; index < fRecPoints->GetEntries(); index++)
142 (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->SetIndexInList(index) ;
143 (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->Print();
148 if(strstr(option,"deb") || strstr(option,"all"))
149 PrintRecPoints(option) ;
151 AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
153 fRecPoints->Delete();
155 if(strstr(option,"tim")){
156 gBenchmark->Stop("EMCALClusterizer");
157 printf("Exec took %f seconds for Clusterizing",
158 gBenchmark->GetCpuTime("EMCALClusterizer"));
162 //____________________________________________________________________________
163 Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
165 // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
167 // = 2 is in different SM; continue searching
168 // In case it is in different SM, but same phi rack, check if neigbours at eta=0
169 // neighbours are defined as digits having at least a common side
170 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
171 // which is compared to a digit (d2) not yet in a cluster
173 static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
174 static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
175 static Int_t rowdiff=0, coldiff=0;
179 fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
180 fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
181 fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
182 fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
184 //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
185 if(nSupMod1 != nSupMod2 )
187 //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
188 Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
189 Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
191 if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
193 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
194 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
195 if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
196 else ieta2+=AliEMCALGeoParams::fgkEMCALCols;
198 shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
200 }//Different SM, same phi
202 rowdiff = TMath::Abs(iphi1 - iphi2);
203 coldiff = TMath::Abs(ieta1 - ieta2) ;
205 // neighbours +-1 in col and row
206 if ( TMath::Abs(coldiff) < 2 && TMath::Abs(rowdiff) < 2)
209 AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
210 d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
216 AliDebug(9, Form("NOT 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));
223 //____________________________________________________________________________
224 void AliEMCALClusterizerNxN::MakeClusters()
226 // Steering method to construct the clusters stored in a list of Reconstructed Points
227 // A cluster is defined as a list of neighbour digits
228 // Mar 03, 2007 by PAI
230 if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
234 // Set up TObjArray with pointers to digits to work on
235 //TObjArray *digitsC = new TObjArray();
237 TIter nextdigit(fDigitsArr);
238 AliEMCALDigit *digit = 0;
239 while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
240 digitsC.AddLast(digit);
243 TIter nextdigitC(&digitsC);
245 AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f)\n",
246 fDigitsArr->GetEntries(),fMinECut));
248 Bool_t bDone = kFALSE;
249 while ( bDone != kTRUE )
251 //first sort the digits:
252 Int_t iMaxEnergyDigit = -1;
253 Float_t dMaxEnergyDigit = -1;
254 AliEMCALDigit *pMaxEnergyDigit = 0;
256 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) )
257 { // scan over the list of digitsC
258 Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
259 //AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
261 //if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold )
262 if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0) // no threshold!
264 if (dEnergyCalibrated > dMaxEnergyDigit)
266 dMaxEnergyDigit = dEnergyCalibrated;
267 iMaxEnergyDigit = digit->GetId();
268 pMaxEnergyDigit = digit;
273 if (iMaxEnergyDigit < 0 || digitsC.GetEntries() <= 0)
279 AliDebug (2, Form("Max digit found: %1.2f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
280 AliDebug(5, Form("Max Digit ENERGY: %1.5f", dMaxEnergyDigit));
282 // keep the candidate digits in a list
283 TList clusterDigitList;
284 clusterDigitList.SetOwner(kFALSE);
285 clusterDigitList.AddLast(pMaxEnergyDigit);
287 Double_t clusterCandidateEnergy = dMaxEnergyDigit;
289 // now loop over the resto of the digits and cluster into NxN cluster
290 // we do not actually cluster yet: we keep them in the list clusterDigitList
292 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) )
293 { // scan over the list of digitsC
294 if (digit == pMaxEnergyDigit) continue;
295 Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
296 AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
297 //if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold )
298 if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0 )
300 Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR?
301 if(TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR?
302 Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
303 if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!!
305 clusterDigitList.AddLast(digit) ;
306 clusterCandidateEnergy += dEnergyCalibrated;
309 }// loop over the next digits
311 // start a cluster here only if a cluster energy is larger than clustering threshold
312 //if (clusterCandidateEnergy > 0.1)
313 AliDebug(5, Form("Clusterization threshold is %f MeV", fECAClusteringThreshold));
314 if (clusterCandidateEnergy > fECAClusteringThreshold)
316 if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
318 AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ;
319 fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
320 recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ;
321 fNumberOfECAClusters++ ;
322 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
324 AliDebug(9, Form("Number of cells per cluster (max is 9!): %d", clusterDigitList.GetEntries()));
325 for (Int_t idig = 0; idig < clusterDigitList.GetEntries(); idig++)
328 digit = (AliEMCALDigit*)clusterDigitList.At(idig);
329 Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
330 AliDebug(5, Form(" Adding digit %d", digit->GetId()));
331 // note: this way the sharing info is lost!
332 recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE) ; //Time or TimeR?
333 digitsC.Remove(digit);
338 // we do not want to start clustering in the same spot!
339 // but in this case we may NOT reuse this seed for another cluster!
340 // need a better bookeeping?
341 digitsC.Remove(pMaxEnergyDigit);
344 AliDebug (2, Form("Number of digits left: %d", digitsC.GetEntries()));
347 //delete digitsC ; //nope we use an object
349 AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
352 //____________________________________________________________________________
353 void AliEMCALClusterizerNxN::MakeUnfolding()
355 // Unfolds clusters using the shape of an ElectroMagnetic shower
356 // Performs unfolding of all clusters
358 if(fNumberOfECAClusters > 0){
360 AliFatal("Did not get geometry from EMCALLoader") ;
361 Int_t nModulesToUnfold = fGeom->GetNCells();
363 Int_t numberofNotUnfolded = fNumberOfECAClusters ;
365 for(index = 0 ; index < numberofNotUnfolded ; index++){
367 AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
371 recPoint->GetGlobalPosition(gpos);
372 fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);
373 if(absId > nModulesToUnfold)
376 Int_t nMultipl = recPoint->GetMultiplicity() ;
377 AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
378 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
379 Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
381 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
382 //UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ;
383 fRecPoints->Remove(recPoint);
384 fRecPoints->Compress() ;
386 fNumberOfECAClusters-- ;
387 numberofNotUnfolded-- ;
390 recPoint->SetNExMax(1) ; //Only one local maximum
394 delete[] maxAtEnergy ;
397 // End of Unfolding of clusters
400 //____________________________________________________________________________
401 Double_t AliEMCALClusterizerNxN::ShowerShape(Double_t x, Double_t y)
403 // Shape of the shower
404 // If you change this function, change also the gradient evaluation in ChiSquare()
406 Double_t r = sqrt(x*x+y*y);
407 Double_t r133 = TMath::Power(r, 1.33) ;
408 Double_t r669 = TMath::Power(r, 6.69) ;
409 Double_t shape = TMath::Exp( -r133 * (1. / (1.57 + 0.0860 * r133) - 0.55 / (1 + 0.000563 * r669) ) ) ;
413 //____________________________________________________________________________
414 void AliEMCALClusterizerNxN::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/,
416 AliEMCALDigit ** /*maxAt*/,
417 Float_t * /*maxAtEnergy*/)
420 // Performs the unfolding of a cluster with nMax overlapping showers
422 AliWarning("Not implemented. To be.");
425 //___________________________________________________________________
426 void AliEMCALClusterizerNxN::SetInputCalibrated(Bool_t val)
429 // input is calibrated - the case when we run already on ESD
431 AliEMCALClusterizerNxN::fgkIsInputCalibrated = val;