]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PHOS/AliPHOSClusterizerv1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Added possibility to accumulate statistics through the number of runs.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSClusterizerv1.cxx
... / ...
CommitLineData
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/* $Id$ */
17
18/* History of cvs commits:
19 *
20 * $Log: AliPHOSClusterizerv1.cxx,v $
21 * Revision 1.118 2007/12/11 21:23:26 kharlov
22 * Added possibility to swith off unfolding
23 *
24 * Revision 1.117 2007/10/18 08:42:05 kharlov
25 * Bad channels cleaned before clusterization
26 *
27 * Revision 1.116 2007/10/01 20:24:08 kharlov
28 * Memory leaks fixed
29 *
30 * Revision 1.115 2007/09/26 14:22:17 cvetan
31 * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
32 *
33 * Revision 1.114 2007/09/06 16:06:44 kharlov
34 * Absence of sorting results in loose of all unfolded clusters
35 *
36 * Revision 1.113 2007/08/28 12:55:07 policheh
37 * Loaders removed from the reconstruction code (C.Cheshkov)
38 *
39 * Revision 1.112 2007/08/22 09:20:50 hristov
40 * Updated QA classes (Yves)
41 *
42 * Revision 1.111 2007/08/08 12:11:28 kharlov
43 * Protection against uninitialized fQADM
44 *
45 * Revision 1.110 2007/08/07 14:16:00 kharlov
46 * Quality assurance added (Yves Schutz)
47 *
48 * Revision 1.109 2007/07/24 17:20:35 policheh
49 * Usage of RecoParam objects instead of hardcoded parameters in reconstruction.
50 * (See $ALICE_ROOT/PHOS/macros/BeamTest2006/RawReconstruction.C).
51 *
52 * Revision 1.108 2007/06/18 07:00:51 kharlov
53 * Bug fix for attempt to use AliPHOSEmcRecPoint after its deletion
54 *
55 * Revision 1.107 2007/05/25 14:12:26 policheh
56 * Local to tracking CS transformation added for CPV rec. points
57 *
58 * Revision 1.106 2007/05/24 13:01:22 policheh
59 * Local to tracking CS transformation invoked for each EMC rec.point
60 *
61 * Revision 1.105 2007/05/02 13:41:22 kharlov
62 * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
63 *
64 * Revision 1.104 2007/04/27 16:55:53 kharlov
65 * Calibration stops if PHOS CDB objects do not exist
66 *
67 * Revision 1.103 2007/04/11 11:55:45 policheh
68 * SetDistancesToBadChannels() added.
69 *
70 * Revision 1.102 2007/03/28 19:18:15 kharlov
71 * RecPoints recalculation in TSM removed
72 *
73 * Revision 1.101 2007/03/06 06:51:27 kharlov
74 * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
75 *
76 * Revision 1.100 2007/01/10 11:07:26 kharlov
77 * Raw digits writing to file (B.Polichtchouk)
78 *
79 * Revision 1.99 2006/11/07 16:49:51 kharlov
80 * Corrections for next event switch in case of raw data (B.Polichtchouk)
81 *
82 * Revision 1.98 2006/10/27 17:14:27 kharlov
83 * Introduce AliDebug and AliLog (B.Polichtchouk)
84 *
85 * Revision 1.97 2006/08/29 11:41:19 kharlov
86 * Missing implementation of ctors and = operator are added
87 *
88 * Revision 1.96 2006/08/25 16:56:30 kharlov
89 * Compliance with Effective C++
90 *
91 * Revision 1.95 2006/08/11 12:36:26 cvetan
92 * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
93 *
94 * Revision 1.94 2006/08/07 12:27:49 hristov
95 * Removing obsolete code which affected the event numbering scheme
96 *
97 * Revision 1.93 2006/08/01 12:20:17 cvetan
98 * 1. Adding a possibility to read and reconstruct an old rcu formatted raw data. This is controlled by an option of AliReconstruction and AliPHOSReconstructor. 2. In case of raw data processing (without galice.root) create the default AliPHOSGeometry object. Most likely this should be moved to the CDB
99 *
100 * Revision 1.92 2006/04/29 20:26:46 hristov
101 * Separate EMC and CPV calibration (Yu.Kharlov)
102 *
103 * Revision 1.91 2006/04/22 10:30:17 hristov
104 * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
105 *
106 * Revision 1.90 2006/04/11 15:22:59 hristov
107 * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
108 *
109 * Revision 1.89 2006/03/13 14:05:42 kharlov
110 * Calibration objects for EMC and CPV
111 *
112 * Revision 1.88 2006/01/11 08:54:52 hristov
113 * Additional protection in case no calibration entry was found
114 *
115 * Revision 1.87 2005/11/22 08:46:43 kharlov
116 * Updated to new CDB (Boris Polichtchouk)
117 *
118 * Revision 1.86 2005/11/14 21:52:43 hristov
119 * Coding conventions
120 *
121 * Revision 1.85 2005/09/27 16:08:08 hristov
122 * New version of CDB storage framework (A.Colla)
123 *
124 * Revision 1.84 2005/09/21 10:02:47 kharlov
125 * Reading calibration from CDB (Boris Polichtchouk)
126 *
127 * Revision 1.82 2005/09/02 15:43:13 kharlov
128 * Add comments in GetCalibrationParameters and Calibrate
129 *
130 * Revision 1.81 2005/09/02 14:32:07 kharlov
131 * Calibration of raw data
132 *
133 * Revision 1.80 2005/08/24 15:31:36 kharlov
134 * Setting raw digits flag
135 *
136 * Revision 1.79 2005/07/25 15:53:53 kharlov
137 * Read raw data
138 *
139 * Revision 1.78 2005/05/28 14:19:04 schutz
140 * Compilation warnings fixed by T.P.
141 *
142 */
143
144//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
145//////////////////////////////////////////////////////////////////////////////
146// Clusterization class. Performs clusterization (collects neighbouring active cells) and
147// unfolds the clusters having several local maxima.
148// Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
149// PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
150// parameters including input digits branch title, thresholds etc.)
151// This TTask is normally called from Reconstructioner, but can as well be used in
152// standalone mode.
153// Use Case:
154// root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1(<pointer_to_phos_geometry_onject>)
155// root [1] cl->Digits2Clusters(digitsTree,clusterTree)
156// //finds RecPoints in the current event
157// root [2] cl->SetDigitsBranch("digits2")
158// //sets another title for Digitis (input) branch
159// root [3] cl->SetRecPointsBranch("recp2")
160// //sets another title four output branches
161// root [4] cl->SetEmcLocalMaxCut(0.03)
162// //set clusterization parameters
163
164// --- ROOT system ---
165
166#include "TMath.h"
167#include "TMinuit.h"
168#include "TTree.h"
169#include "TBenchmark.h"
170#include "TClonesArray.h"
171
172// --- Standard library ---
173
174// --- AliRoot header files ---
175#include "AliLog.h"
176#include "AliConfig.h"
177#include "AliPHOSGeometry.h"
178#include "AliPHOSClusterizerv1.h"
179#include "AliPHOSEmcRecPoint.h"
180#include "AliPHOSCpvRecPoint.h"
181#include "AliPHOSDigit.h"
182#include "AliPHOSDigitizer.h"
183#include "AliCDBManager.h"
184#include "AliCDBStorage.h"
185#include "AliCDBEntry.h"
186#include "AliPHOSRecoParam.h"
187#include "AliPHOSReconstructor.h"
188#include "AliPHOSCalibData.h"
189
190ClassImp(AliPHOSClusterizerv1)
191
192//____________________________________________________________________________
193AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
194 AliPHOSClusterizer(),
195 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
196 fWrite(0),
197 fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
198 fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
199 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
200 fW0CPV(0), fEmcTimeGate(0), fEcoreRadius(0)
201{
202 // default ctor (to be used mainly by Streamer)
203
204 fDefaultInit = kTRUE ;
205}
206
207//____________________________________________________________________________
208AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
209 AliPHOSClusterizer(geom),
210 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
211 fWrite(0),
212 fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
213 fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
214 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
215 fW0CPV(0), fEmcTimeGate(0), fEcoreRadius(0)
216{
217 // ctor with the indication of the file where header Tree and digits Tree are stored
218
219 Init() ;
220 fDefaultInit = kFALSE ;
221}
222
223//____________________________________________________________________________
224 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
225{
226 // dtor
227
228}
229//____________________________________________________________________________
230void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
231{
232 // Steering method to perform clusterization for one event
233 // The input is the tree with digits.
234 // The output is the tree with clusters.
235
236 if(strstr(option,"tim"))
237 gBenchmark->Start("PHOSClusterizer");
238
239 if(strstr(option,"print")) {
240 Print() ;
241 return ;
242 }
243
244 MakeClusters() ;
245
246 AliDebug(2,Form(" ---- Printing clusters (%d)\n",
247 fEMCRecPoints->GetEntries()));
248 if(AliLog::GetGlobalDebugLevel()>1)
249 fEMCRecPoints->Print();
250
251 if(fToUnfold)
252 MakeUnfolding();
253
254 WriteRecPoints();
255
256 if(strstr(option,"deb"))
257 PrintRecPoints(option) ;
258
259 if(strstr(option,"tim")){
260 gBenchmark->Stop("PHOSClusterizer");
261 AliInfo(Form("took %f seconds for Clusterizing\n",
262 gBenchmark->GetCpuTime("PHOSClusterizer")));
263 }
264 fEMCRecPoints->Delete();
265 fCPVRecPoints->Delete();
266}
267
268//____________________________________________________________________________
269Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
270 Int_t nPar, Float_t * fitparameters) const
271{
272 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
273 // The initial values for fitting procedure are set equal to the positions of local maxima.
274 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
275
276
277 if(!gMinuit) //it was deleted by someone else
278 gMinuit = new TMinuit(100) ;
279 gMinuit->mncler(); // Reset Minuit's list of paramters
280 gMinuit->SetPrintLevel(-1) ; // No Printout
281 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
282 // To set the address of the minimization function
283
284 TList * toMinuit = new TList();
285 toMinuit->AddAt(emcRP,0) ;
286 toMinuit->AddAt(fDigitsArr,1) ;
287 toMinuit->AddAt(fGeom,2) ;
288
289 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
290
291 // filling initial values for fit parameters
292 AliPHOSDigit * digit ;
293
294 Int_t ierflg = 0;
295 Int_t index = 0 ;
296 Int_t nDigits = (Int_t) nPar / 3 ;
297
298 Int_t iDigit ;
299
300 for(iDigit = 0; iDigit < nDigits; iDigit++){
301 digit = maxAt[iDigit];
302
303 Int_t relid[4] ;
304 Float_t x = 0.;
305 Float_t z = 0.;
306 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
307 fGeom->RelPosInModule(relid, x, z) ;
308
309 Float_t energy = maxAtEnergy[iDigit] ;
310
311 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
312 index++ ;
313 if(ierflg != 0){
314 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
315 return kFALSE;
316 }
317 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
318 index++ ;
319 if(ierflg != 0){
320 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
321 return kFALSE;
322 }
323 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
324 index++ ;
325 if(ierflg != 0){
326 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
327 return kFALSE;
328 }
329 }
330
331 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
332 // depends on it.
333 Double_t p1 = 1.0 ;
334 Double_t p2 = 0.0 ;
335
336 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
337 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
338 gMinuit->SetMaxIterations(5);
339 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
340
341 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
342
343 if(ierflg == 4){ // Minimum not found
344 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
345 return kFALSE ;
346 }
347 for(index = 0; index < nPar; index++){
348 Double_t err ;
349 Double_t val ;
350 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
351 fitparameters[index] = val ;
352 }
353
354 delete toMinuit ;
355 return kTRUE;
356
357}
358
359
360//____________________________________________________________________________
361void AliPHOSClusterizerv1::Init()
362{
363 // Make all memory allocations which can not be done in default constructor.
364 // Attach the Clusterizer task to the list of PHOS tasks
365
366 fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
367
368 if(!gMinuit)
369 gMinuit = new TMinuit(100);
370
371 if (!fgCalibData)
372 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
373 if (fgCalibData->GetCalibDataEmc() == 0)
374 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
375 if (fgCalibData->GetCalibDataCpv() == 0)
376 AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
377
378}
379
380//____________________________________________________________________________
381void AliPHOSClusterizerv1::InitParameters()
382{
383
384 fNumberOfCpvClusters = 0 ;
385 fNumberOfEmcClusters = 0 ;
386
387 const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
388 if(!recoParam) AliFatal("Reconstruction parameters are not set!");
389
390 recoParam->Print();
391
392 fEmcClusteringThreshold = recoParam->GetEMCClusteringThreshold();
393 fCpvClusteringThreshold = recoParam->GetCPVClusteringThreshold();
394
395 fEmcLocMaxCut = recoParam->GetEMCLocalMaxCut();
396 fCpvLocMaxCut = recoParam->GetCPVLocalMaxCut();
397
398 fW0 = recoParam->GetEMCLogWeight();
399 fW0CPV = recoParam->GetCPVLogWeight();
400
401 fEmcTimeGate = 1.e-6 ; //10 sample steps
402 fEcoreRadius = recoParam->GetEMCEcoreRadius();
403
404 fToUnfold = recoParam->EMCToUnfold() ;
405
406 fWrite = kTRUE ;
407}
408
409//____________________________________________________________________________
410Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
411{
412 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
413 // = 1 are neighbour
414 // = 2 are not neighbour but do not continue searching
415 // =-1 are not neighbour, continue searching, but do not look before d2 next time
416 // neighbours are defined as digits having at least a common vertex
417 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
418 // which is compared to a digit (d2) not yet in a cluster
419
420 Int_t relid1[4] ;
421 fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
422
423 Int_t relid2[4] ;
424 fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
425
426 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
427 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
428 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
429
430 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
431 // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
432 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
433 return 1 ;
434 }
435 else {
436 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
437 return 2; // Difference in row numbers is too large to look further
438 }
439 return 0 ;
440
441 }
442 else {
443 if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
444 return -1 ;
445 if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
446 return -1 ;
447
448 return 2 ;
449
450 }
451
452 return 0 ;
453}
454//____________________________________________________________________________
455Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
456{
457 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
458
459 Bool_t rv = kFALSE ;
460
461 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
462
463 if(digit->GetId() <= nEMC ) rv = kTRUE;
464
465 return rv ;
466}
467
468//____________________________________________________________________________
469Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
470{
471 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
472
473 Bool_t rv = kFALSE ;
474
475 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
476
477 if(digit->GetId() > nEMC ) rv = kTRUE;
478
479 return rv ;
480}
481
482//____________________________________________________________________________
483void AliPHOSClusterizerv1::WriteRecPoints()
484{
485
486 // Creates new branches with given title
487 // fills and writes into TreeR.
488
489 Int_t index ;
490 //Evaluate position, dispersion and other RecPoint properties..
491 Int_t nEmc = fEMCRecPoints->GetEntriesFast();
492 Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
493 TVector3 fakeVtx(0.,0.,0.) ;
494 for(index = 0; index < nEmc; index++){
495 AliPHOSEmcRecPoint * rp =
496 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
497 rp->Purify(emcMinE) ;
498 if(rp->GetMultiplicity()==0){
499 fEMCRecPoints->RemoveAt(index) ;
500 delete rp ;
501 continue;
502 }
503
504 // No vertex is available now, calculate corrections in PID
505 rp->EvalAll(fDigitsArr) ;
506 rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
507 rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
508 rp->EvalLocal2TrackingCSTransform();
509 }
510 fEMCRecPoints->Compress() ;
511 fEMCRecPoints->Sort() ;
512 // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
513 for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
514 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
515 }
516
517 //For each rec.point set the distance to the nearest bad crystal (BVP)
518 SetDistancesToBadChannels();
519
520 //Now the same for CPV
521 for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
522 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
523 rp->EvalAll(fDigitsArr) ;
524 rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
525 rp->EvalLocal2TrackingCSTransform();
526 }
527 fCPVRecPoints->Sort() ;
528
529 for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
530 dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
531
532 fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
533
534 if(fWrite){ //We write TreeR
535 fTreeR->Fill();
536 }
537}
538
539//____________________________________________________________________________
540void AliPHOSClusterizerv1::MakeClusters()
541{
542 // Steering method to construct the clusters stored in a list of Reconstructed Points
543 // A cluster is defined as a list of neighbour digits
544
545 fNumberOfCpvClusters = 0 ;
546 fNumberOfEmcClusters = 0 ;
547
548 //Mark all digits as unused yet
549 const Int_t maxNDigits = 3584; // There is no clusters larger than PHOS module ;)
550 Int_t nDigits=fDigitsArr->GetEntriesFast() ;
551
552 for(Int_t i=0; i<nDigits; i++){
553 fDigitsUsed[i]=0 ;
554 }
555 Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
556 //e.g. first digit in this module, first CPV digit etc.
557 AliPHOSDigit * digit ;
558 TArrayI clusterdigitslist(maxNDigits) ;
559 AliPHOSRecPoint * clu = 0 ;
560 for(Int_t i=0; i<nDigits; i++){
561 if(fDigitsUsed[i])
562 continue ;
563
564 digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
565
566 clu=0 ;
567
568 Int_t index ;
569
570 //is this digit so energetic that start cluster?
571 if (( IsInEmc(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
572 ( IsInCpv(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fCpvClusteringThreshold ) ) {
573 Int_t iDigitInCluster = 0 ;
574 if ( IsInEmc(digit) ) {
575 // start a new EMC RecPoint
576 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
577 fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
578
579 fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
580 clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
581 fNumberOfEmcClusters++ ;
582 clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),CalibrateT(digit->GetTime(),digit->GetId())) ;
583 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
584 iDigitInCluster++ ;
585 fDigitsUsed[i]=kTRUE ;
586 } else {
587 // start a new CPV cluster
588 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
589 fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
590
591 fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
592 clu = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
593 fNumberOfCpvClusters++ ;
594 clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),0.) ; // no timing information in CPV
595 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
596 iDigitInCluster++ ;
597 fDigitsUsed[i]=kTRUE ;
598 } // else
599
600 //Now scan remaining digits in list to find neigbours of our seed
601
602 AliPHOSDigit * digitN ;
603 index = 0 ;
604 while (index < iDigitInCluster){ // scan over digits already in cluster
605 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
606 index++ ;
607 for(Int_t j=iFirst; j<nDigits; j++){
608 if (iDigitInCluster >= maxNDigits) {
609 AliError(Form("The number of digits in cluster is more than %d, skip the rest of event",
610 maxNDigits));
611 return;
612 }
613 if(fDigitsUsed[j])
614 continue ; //look through remaining digits
615 digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
616 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
617 switch (ineb ) {
618 case -1: //too early (e.g. previous module), do not look before j at subsequent passes
619 iFirst=j ;
620 break ;
621 case 0 : // not a neighbour
622 break ;
623 case 1 : // are neighbours
624 clu->AddDigit(*digitN, Calibrate(digitN->GetEnergy(),digit->GetId()),CalibrateT(digitN->GetTime(),digit->GetId())) ;
625 clusterdigitslist[iDigitInCluster] = j ;
626 iDigitInCluster++ ;
627 fDigitsUsed[j]=kTRUE ;
628 break ;
629 case 2 : // too far from each other
630 goto endOfLoop;
631 } // switch
632
633 }
634
635 endOfLoop: ; //scanned all possible neighbours for this digit
636
637 } // loop over cluster
638 } // energy theshold
639 }
640
641}
642
643//____________________________________________________________________________
644void AliPHOSClusterizerv1::MakeUnfolding()
645{
646 // Unfolds clusters using the shape of an ElectroMagnetic shower
647 // Performs unfolding of all EMC/CPV clusters
648
649 // Unfold first EMC clusters
650 if(fNumberOfEmcClusters > 0){
651
652 Int_t nModulesToUnfold = fGeom->GetNModules() ;
653
654 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
655 Int_t index ;
656 for(index = 0 ; index < numberofNotUnfolded ; index++){
657
658 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
659 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
660 break ;
661
662 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
663 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
664 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
665 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
666
667 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
668 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
669
670 fEMCRecPoints->Remove(emcRecPoint);
671 fEMCRecPoints->Compress() ;
672 index-- ;
673 fNumberOfEmcClusters -- ;
674 numberofNotUnfolded-- ;
675 }
676 else{
677 emcRecPoint->SetNExMax(1) ; //Only one local maximum
678 }
679
680 delete[] maxAt ;
681 delete[] maxAtEnergy ;
682 }
683 }
684 // Unfolding of EMC clusters finished
685
686
687 // Unfold now CPV clusters
688 if(fNumberOfCpvClusters > 0){
689
690 Int_t nModulesToUnfold = fGeom->GetNModules() ;
691
692 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
693 Int_t index ;
694 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
695
696 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
697
698 if(recPoint->GetPHOSMod()> nModulesToUnfold)
699 break ;
700
701 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
702
703 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
704 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
705 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
706 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
707
708 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
709 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
710 fCPVRecPoints->Remove(emcRecPoint);
711 fCPVRecPoints->Compress() ;
712 index-- ;
713 numberofCpvNotUnfolded-- ;
714 fNumberOfCpvClusters-- ;
715 }
716
717 delete[] maxAt ;
718 delete[] maxAtEnergy ;
719 }
720 }
721 //Unfolding of Cpv clusters finished
722
723}
724
725//____________________________________________________________________________
726Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
727{
728 // Shape of the shower (see PHOS TDR)
729 // If you change this function, change also the gradient evaluation in ChiSquare()
730
731 //for the moment we neglect dependence on the incident angle.
732
733 Double_t r2 = x*x + z*z ;
734 Double_t r4 = r2*r2 ;
735 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
736 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
737 return shape ;
738}
739
740//____________________________________________________________________________
741void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
742 Int_t nMax,
743 AliPHOSDigit ** maxAt,
744 Float_t * maxAtEnergy)
745{
746 // Performs the unfolding of a cluster with nMax overlapping showers
747
748 Int_t nPar = 3 * nMax ;
749 Float_t * fitparameters = new Float_t[nPar] ;
750
751 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
752
753 if( !rv ) {
754 // Fit failed, return and remove cluster
755 iniEmc->SetNExMax(-1) ;
756 delete[] fitparameters ;
757 return ;
758 }
759
760 // create ufolded rec points and fill them with new energy lists
761 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
762 // and later correct this number in acordance with actual energy deposition
763
764 Int_t nDigits = iniEmc->GetMultiplicity() ;
765 Float_t * efit = new Float_t[nDigits] ;
766 Float_t xDigit=0.,zDigit=0. ;
767 Float_t xpar=0.,zpar=0.,epar=0. ;
768 Int_t relid[4] ;
769 AliPHOSDigit * digit = 0 ;
770 Int_t * emcDigits = iniEmc->GetDigitsList() ;
771
772 TVector3 vIncid ;
773
774 Int_t iparam ;
775 Int_t iDigit ;
776 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
777 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
778 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
779 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
780 efit[iDigit] = 0;
781
782 iparam = 0 ;
783 while(iparam < nPar ){
784 xpar = fitparameters[iparam] ;
785 zpar = fitparameters[iparam+1] ;
786 epar = fitparameters[iparam+2] ;
787 iparam += 3 ;
788// fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
789// efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
790 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
791 }
792 }
793
794 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
795 // so that energy deposited in each cell is distributed betwin new clusters proportionally
796 // to its contribution to efit
797
798 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
799 Float_t ratio ;
800
801 iparam = 0 ;
802 while(iparam < nPar ){
803 xpar = fitparameters[iparam] ;
804 zpar = fitparameters[iparam+1] ;
805 epar = fitparameters[iparam+2] ;
806 iparam += 3 ;
807// fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
808
809 AliPHOSEmcRecPoint * emcRP = 0 ;
810
811 if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
812
813 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
814 fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
815
816 (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
817 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
818 fNumberOfEmcClusters++ ;
819 emcRP->SetNExMax((Int_t)nPar/3) ;
820 }
821 else{//create new entries in fCPVRecPoints
822 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
823 fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
824
825 (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
826 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
827 fNumberOfCpvClusters++ ;
828 }
829
830 Float_t eDigit ;
831 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
832 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
833 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
834 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
835// ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
836 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
837 eDigit = emcEnergies[iDigit] * ratio ;
838 emcRP->AddDigit( *digit, eDigit,CalibrateT(digit->GetTime(),digit->GetId()) ) ;
839 }
840 }
841
842 delete[] fitparameters ;
843 delete[] efit ;
844
845}
846
847//_____________________________________________________________________________
848void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
849{
850 // Calculates the Chi square for the cluster unfolding minimization
851 // Number of parameters, Gradient, Chi squared, parameters, what to do
852
853 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
854
855 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
856 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
857 // A bit buggy way to get an access to the geometry
858 // To be revised!
859 AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
860
861// TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
862
863 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
864
865 Int_t * emcDigits = emcRP->GetDigitsList() ;
866
867 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
868
869 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
870
871// TVector3 vInc ;
872 fret = 0. ;
873 Int_t iparam ;
874
875 if(iflag == 2)
876 for(iparam = 0 ; iparam < nPar ; iparam++)
877 Grad[iparam] = 0 ; // Will evaluate gradient
878
879 Double_t efit ;
880
881 AliPHOSDigit * digit ;
882 Int_t iDigit ;
883
884 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
885
886 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
887
888 Int_t relid[4] ;
889 Float_t xDigit ;
890 Float_t zDigit ;
891
892 geom->AbsToRelNumbering(digit->GetId(), relid) ;
893
894 geom->RelPosInModule(relid, xDigit, zDigit) ;
895
896 if(iflag == 2){ // calculate gradient
897 Int_t iParam = 0 ;
898 efit = 0 ;
899 while(iParam < nPar ){
900 Double_t dx = (xDigit - x[iParam]) ;
901 iParam++ ;
902 Double_t dz = (zDigit - x[iParam]) ;
903 iParam++ ;
904// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
905// efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
906 efit += x[iParam] * ShowerShape(dx,dz) ;
907 iParam++ ;
908 }
909 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
910 iParam = 0 ;
911 while(iParam < nPar ){
912 Double_t xpar = x[iParam] ;
913 Double_t zpar = x[iParam+1] ;
914 Double_t epar = x[iParam+2] ;
915// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
916 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
917// Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
918 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
919//DP: No incident angle dependence in gradient yet!!!!!!
920 Double_t r4 = dr*dr*dr*dr ;
921 Double_t r295 = TMath::Power(dr,2.95) ;
922 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
923 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
924
925 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
926 iParam++ ;
927 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
928 iParam++ ;
929 Grad[iParam] += shape ; // Derivative over energy
930 iParam++ ;
931 }
932 }
933 efit = 0;
934 iparam = 0 ;
935
936 while(iparam < nPar ){
937 Double_t xpar = x[iparam] ;
938 Double_t zpar = x[iparam+1] ;
939 Double_t epar = x[iparam+2] ;
940 iparam += 3 ;
941// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
942// efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
943 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
944 }
945
946 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
947 // Here we assume, that sigma = sqrt(E)
948 }
949
950}
951
952//____________________________________________________________________________
953void AliPHOSClusterizerv1::Print(const Option_t *)const
954{
955 // Print clusterizer parameters
956
957 TString message ;
958 TString taskName(GetName()) ;
959 taskName.ReplaceAll(Version(), "") ;
960
961 if( strcmp(GetName(), "") !=0 ) {
962 // Print parameters
963 message = "\n--------------- %s %s -----------\n" ;
964 message += "Clusterizing digits from the file: %s\n" ;
965 message += " Branch: %s\n" ;
966 message += " EMC Clustering threshold = %f\n" ;
967 message += " EMC Local Maximum cut = %f\n" ;
968 message += " EMC Logarothmic weight = %f\n" ;
969 message += " CPV Clustering threshold = %f\n" ;
970 message += " CPV Local Maximum cut = %f\n" ;
971 message += " CPV Logarothmic weight = %f\n" ;
972 if(fToUnfold)
973 message += " Unfolding on\n" ;
974 else
975 message += " Unfolding off\n" ;
976
977 message += "------------------------------------------------------------------" ;
978 }
979 else
980 message = " AliPHOSClusterizerv1 not initialized " ;
981
982 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
983 taskName.Data(),
984 GetTitle(),
985 taskName.Data(),
986 GetName(),
987 fEmcClusteringThreshold,
988 fEmcLocMaxCut,
989 fW0,
990 fCpvClusteringThreshold,
991 fCpvLocMaxCut,
992 fW0CPV )) ;
993}
994//____________________________________________________________________________
995void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
996{
997 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
998
999 AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
1000 fEMCRecPoints->GetEntriesFast(),
1001 fCPVRecPoints->GetEntriesFast() )) ;
1002
1003 if(strstr(option,"all")) {
1004 printf("\n EMC clusters \n") ;
1005 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1006 Int_t index ;
1007 for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1008 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1009 TVector3 locpos;
1010 rp->GetLocalPosition(locpos);
1011 Float_t lambda[2];
1012 rp->GetElipsAxis(lambda);
1013 Int_t * primaries;
1014 Int_t nprimaries;
1015 primaries = rp->GetPrimaries(nprimaries);
1016 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1017 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1018 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1019
1020 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1021 printf("%d ", primaries[iprimary] ) ;
1022 }
1023 printf("\n") ;
1024 }
1025
1026 //Now plot CPV recPoints
1027 printf("\n CPV clusters \n") ;
1028 printf("Index Ene(MeV) Module X Y Z \n") ;
1029 for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1030 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1031
1032 TVector3 locpos;
1033 rp->GetLocalPosition(locpos);
1034
1035 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1036 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1037 locpos.X(), locpos.Y(), locpos.Z()) ;
1038 }
1039 }
1040}
1041
1042
1043//____________________________________________________________________________
1044void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1045{
1046 //For each EMC rec. point set the distance to the nearest bad crystal.
1047 //Author: Boris Polichtchouk
1048
1049 if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1050
1051 Int_t badIds[8000];
1052 fgCalibData->EmcBadChannelIds(badIds);
1053
1054 AliPHOSEmcRecPoint* rp;
1055
1056 TMatrixF gmat;
1057 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1058 TVector3 gposBadChannel; // global position of bad crystal
1059 TVector3 dR;
1060
1061 Float_t dist,minDist;
1062 Int_t relid[4]={0,0,0,0} ;
1063 TVector3 lpos ;
1064 for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1065 rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1066 //evaluate distance to border
1067 relid[0]=rp->GetPHOSMod() ;
1068 relid[2]=1 ;
1069 relid[3]=1 ;
1070 Float_t xcorner,zcorner;
1071 fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1072 rp->GetLocalPosition(lpos) ;
1073 minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
1074 for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1075 fGeom->AbsToRelNumbering(badIds[iBad],relid) ;
1076 if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1077 continue ; //bad channels can be in the module which does not exist in simulations.
1078 rp->GetGlobalPosition(gposRecPoint,gmat);
1079 fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1080 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1081 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1082 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1083 dR = gposBadChannel-gposRecPoint;
1084 dist = dR.Mag();
1085 if(dist<minDist) minDist = dist;
1086 }
1087
1088 rp->SetDistanceToBadCrystal(minDist);
1089 }
1090
1091}
1092//==================================================================================
1093Float_t AliPHOSClusterizerv1::Calibrate(Float_t amp, Int_t absId) const{
1094 // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB
1095
1096 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1097
1098 //Determine rel.position of the cell absolute ID
1099 Int_t relId[4];
1100 geom->AbsToRelNumbering(absId,relId);
1101 Int_t module=relId[0];
1102 Int_t row =relId[2];
1103 Int_t column=relId[3];
1104 if(relId[1]){ //CPV
1105 Float_t calibration = fgCalibData->GetADCchannelCpv(module,column,row);
1106 return amp*calibration ;
1107 }
1108 else{ //EMC
1109 Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row);
1110 return amp*calibration ;
1111 }
1112}
1113//==================================================================================
1114Float_t AliPHOSClusterizerv1::CalibrateT(Float_t time, Int_t absId)const{
1115 // Calibrate time in EMC digit
1116
1117 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1118
1119 //Determine rel.position of the cell absolute ID
1120 Int_t relId[4];
1121 geom->AbsToRelNumbering(absId,relId);
1122 Int_t module=relId[0];
1123 Int_t row =relId[2];
1124 Int_t column=relId[3];
1125 if(relId[1]){ //CPV
1126 return 0. ;
1127 }
1128 else{ //EMC
1129 time += fgCalibData->GetTimeShiftEmc(module,column,row);
1130 return time ;
1131 }
1132}
1133