]> 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
TOF pp spectra task added to PWGLFspectra lib
[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),
201 fTimeGateLowAmp(0.), fTimeGateLow(0.), fTimeGateHigh(0.),
202 fEcoreRadius(0)
203{
204 // default ctor (to be used mainly by Streamer)
205
206 fDefaultInit = kTRUE ;
207
208 for(Int_t i=0; i<53760; i++){
209 fDigitsUsed[i]=0 ;
210 }
211}
212
213//____________________________________________________________________________
214AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
215 AliPHOSClusterizer(geom),
216 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
217 fWrite(0),
218 fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
219 fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
220 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
221 fW0CPV(0),
222 fTimeGateLowAmp(0.), fTimeGateLow(0.), fTimeGateHigh(0.),
223 fEcoreRadius(0)
224{
225 // ctor with the indication of the file where header Tree and digits Tree are stored
226
227 for(Int_t i=0; i<53760; i++){
228 fDigitsUsed[i]=0 ;
229 }
230
231 Init() ;
232 fDefaultInit = kFALSE ;
233}
234
235//____________________________________________________________________________
236 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
237{
238 // dtor
239
240}
241//____________________________________________________________________________
242void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
243{
244 // Steering method to perform clusterization for one event
245 // The input is the tree with digits.
246 // The output is the tree with clusters.
247
248 if(strstr(option,"tim"))
249 gBenchmark->Start("PHOSClusterizer");
250
251 if(strstr(option,"print")) {
252 Print() ;
253 return ;
254 }
255
256 MakeClusters() ;
257
258 AliDebug(2,Form(" ---- Printing clusters (%d)\n",
259 fEMCRecPoints->GetEntries()));
260 if(AliLog::GetGlobalDebugLevel()>1)
261 fEMCRecPoints->Print();
262
263 if(fToUnfold)
264 MakeUnfolding();
265
266 WriteRecPoints();
267
268 if(strstr(option,"deb"))
269 PrintRecPoints(option) ;
270
271 if(strstr(option,"tim")){
272 gBenchmark->Stop("PHOSClusterizer");
273 AliInfo(Form("took %f seconds for Clusterizing\n",
274 gBenchmark->GetCpuTime("PHOSClusterizer")));
275 }
276 fEMCRecPoints->Delete();
277 fCPVRecPoints->Delete();
278}
279
280//____________________________________________________________________________
281Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
282 Int_t nPar, Float_t * fitparameters) const
283{
284 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
285 // The initial values for fitting procedure are set equal to the positions of local maxima.
286 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
287
288
289 if(!gMinuit) //it was deleted by someone else
290 gMinuit = new TMinuit(100) ;
291 gMinuit->mncler(); // Reset Minuit's list of paramters
292 gMinuit->SetPrintLevel(-1) ; // No Printout
293 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
294 // To set the address of the minimization function
295
296 TList * toMinuit = new TList();
297 toMinuit->AddAt(emcRP,0) ;
298 toMinuit->AddAt(fDigitsArr,1) ;
299 toMinuit->AddAt(fGeom,2) ;
300
301 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
302
303 // filling initial values for fit parameters
304 AliPHOSDigit * digit ;
305
306 Int_t ierflg = 0;
307 Int_t index = 0 ;
308 Int_t nDigits = (Int_t) nPar / 3 ;
309
310 Int_t iDigit ;
311
312 for(iDigit = 0; iDigit < nDigits; iDigit++){
313 digit = maxAt[iDigit];
314
315 Int_t relid[4] ;
316 Float_t x = 0.;
317 Float_t z = 0.;
318 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
319 fGeom->RelPosInModule(relid, x, z) ;
320
321 Float_t energy = maxAtEnergy[iDigit] ;
322
323 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
324 index++ ;
325 if(ierflg != 0){
326 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
327 return kFALSE;
328 }
329 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
330 index++ ;
331 if(ierflg != 0){
332 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
333 return kFALSE;
334 }
335 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
336 index++ ;
337 if(ierflg != 0){
338 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
339 return kFALSE;
340 }
341 }
342
343 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
344 // depends on it.
345 Double_t p1 = 1.0 ;
346 Double_t p2 = 0.0 ;
347
348 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
349 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
350 gMinuit->SetMaxIterations(5);
351 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
352
353 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
354
355 if(ierflg == 4){ // Minimum not found
356 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
357 return kFALSE ;
358 }
359 for(index = 0; index < nPar; index++){
360 Double_t err ;
361 Double_t val ;
362 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
363 fitparameters[index] = val ;
364 }
365
366 delete toMinuit ;
367 return kTRUE;
368
369}
370
371
372//____________________________________________________________________________
373void AliPHOSClusterizerv1::Init()
374{
375 // Make all memory allocations which can not be done in default constructor.
376 // Attach the Clusterizer task to the list of PHOS tasks
377
378 fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
379
380 if(!gMinuit)
381 gMinuit = new TMinuit(100);
382
383 if (!fgCalibData)
384 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
385 if (fgCalibData->GetCalibDataEmc() == 0)
386 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
387 if (fgCalibData->GetCalibDataCpv() == 0)
388 AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
389
390}
391
392//____________________________________________________________________________
393void AliPHOSClusterizerv1::InitParameters()
394{
395
396 fNumberOfCpvClusters = 0 ;
397 fNumberOfEmcClusters = 0 ;
398
399 const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
400 if(!recoParam) AliFatal("Reconstruction parameters are not set!");
401
402 recoParam->Print();
403
404 fEmcClusteringThreshold = recoParam->GetEMCClusteringThreshold();
405 fCpvClusteringThreshold = recoParam->GetCPVClusteringThreshold();
406
407 fEmcLocMaxCut = recoParam->GetEMCLocalMaxCut();
408 fCpvLocMaxCut = recoParam->GetCPVLocalMaxCut();
409
410 fW0 = recoParam->GetEMCLogWeight();
411 fW0CPV = recoParam->GetCPVLogWeight();
412
413 fTimeGateLowAmp = recoParam->GetTimeGateAmpThresh() ;
414 fTimeGateLow = recoParam->GetTimeGateLow() ;
415 fTimeGateHigh = recoParam->GetTimeGateHigh() ;
416
417 fEcoreRadius = recoParam->GetEMCEcoreRadius();
418
419 fToUnfold = recoParam->EMCToUnfold() ;
420
421 fWrite = kTRUE ;
422}
423
424//____________________________________________________________________________
425Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
426{
427 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
428 // = 1 are neighbour
429 // = 2 are not neighbour but do not continue searching
430 // =-1 are not neighbour, continue searching, but do not look before d2 next time
431 // neighbours are defined as digits having at least a common vertex
432 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
433 // which is compared to a digit (d2) not yet in a cluster
434
435 Int_t relid1[4] ;
436 fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
437
438 Int_t relid2[4] ;
439 fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
440
441 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
442 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
443 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
444
445 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
446 // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
447 if((relid1[1] != 0) || CheckTimeGate(d1->GetTime(),d1->GetEnergy(),d2->GetTime(),d2->GetEnergy()))
448 return 1 ;
449 }
450 else {
451 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
452 return 2; // Difference in row numbers is too large to look further
453 }
454 return 0 ;
455
456 }
457 else {
458 if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
459 return -1 ;
460 if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
461 return -1 ;
462
463 return 2 ;
464
465 }
466
467 return 0 ;
468}
469//____________________________________________________________________________
470Bool_t AliPHOSClusterizerv1::CheckTimeGate(Float_t t1, Float_t amp1, Float_t t2, Float_t amp2)const{
471 //Check if two cells have reasonable time difference
472 //Note that at low amplitude time is defined up to 1 tick == 100 ns.
473 if(amp1<fTimeGateLowAmp || amp2<fTimeGateLowAmp){
474 return (TMath::Abs(t1 - t2 ) < fTimeGateLow) ;
475 }
476 else{ //Time should be measured with good accuracy
477 return (TMath::Abs(t1 - t2 ) < fTimeGateHigh) ;
478 }
479
480}
481//____________________________________________________________________________
482Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
483{
484 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
485
486 Bool_t rv = kFALSE ;
487
488 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
489
490 if(digit->GetId() <= nEMC ) rv = kTRUE;
491
492 return rv ;
493}
494
495//____________________________________________________________________________
496Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
497{
498 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
499
500 Bool_t rv = kFALSE ;
501
502 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
503
504 if(digit->GetId() > nEMC ) rv = kTRUE;
505
506 return rv ;
507}
508
509//____________________________________________________________________________
510void AliPHOSClusterizerv1::WriteRecPoints()
511{
512
513 // Creates new branches with given title
514 // fills and writes into TreeR.
515
516 Int_t index ;
517 //Evaluate position, dispersion and other RecPoint properties..
518 Int_t nEmc = fEMCRecPoints->GetEntriesFast();
519 Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
520 TVector3 fakeVtx(0.,0.,0.) ;
521 for(index = 0; index < nEmc; index++){
522 AliPHOSEmcRecPoint * rp =
523 static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
524 rp->Purify(emcMinE) ;
525 if(rp->GetMultiplicity()==0){
526 fEMCRecPoints->RemoveAt(index) ;
527 delete rp ;
528 continue;
529 }
530
531 // No vertex is available now, calculate corrections in PID
532 rp->EvalAll(fDigitsArr) ;
533 rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
534 rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
535 rp->EvalLocal2TrackingCSTransform();
536 }
537 fEMCRecPoints->Compress() ;
538 fEMCRecPoints->Sort() ;
539 // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
540 for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
541 static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
542 }
543
544 //For each rec.point set the distance to the nearest bad crystal (BVP)
545 SetDistancesToBadChannels();
546
547 //Now the same for CPV
548 for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
549 AliPHOSCpvRecPoint * rp = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
550 rp->EvalAll(fDigitsArr) ;
551 rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
552 rp->EvalLocal2TrackingCSTransform();
553 }
554 fCPVRecPoints->Sort() ;
555
556 for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
557 static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
558
559 fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
560
561 if(fWrite){ //We write TreeR
562 fTreeR->Fill();
563 }
564}
565
566//____________________________________________________________________________
567void AliPHOSClusterizerv1::MakeClusters()
568{
569 // Steering method to construct the clusters stored in a list of Reconstructed Points
570 // A cluster is defined as a list of neighbour digits
571
572 fNumberOfCpvClusters = 0 ;
573 fNumberOfEmcClusters = 0 ;
574
575 //Mark all digits as unused yet
576 const Int_t maxNDigits = 3584; // There is no clusters larger than PHOS module ;)
577 Int_t nDigits=fDigitsArr->GetEntriesFast() ;
578
579 for(Int_t i=0; i<nDigits; i++){
580 fDigitsUsed[i]=0 ;
581 }
582 Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
583 //e.g. first digit in this module, first CPV digit etc.
584 AliPHOSDigit * digit ;
585 TArrayI clusterdigitslist(maxNDigits) ;
586 AliPHOSRecPoint * clu = 0 ;
587 for(Int_t i=0; i<nDigits; i++){
588 if(fDigitsUsed[i])
589 continue ;
590
591 digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
592
593 clu=0 ;
594
595 Int_t index ;
596
597 //is this digit so energetic that start cluster?
598 if (( IsInEmc(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
599 ( IsInCpv(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fCpvClusteringThreshold ) ) {
600 Int_t iDigitInCluster = 0 ;
601 if ( IsInEmc(digit) ) {
602 // start a new EMC RecPoint
603 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
604 fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
605
606 fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
607 clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
608 fNumberOfEmcClusters++ ;
609 clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),CalibrateT(digit->GetTime(),digit->GetId())) ;
610 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
611 iDigitInCluster++ ;
612 fDigitsUsed[i]=kTRUE ;
613 } else {
614 // start a new CPV cluster
615 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
616 fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
617
618 fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
619 clu = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
620 fNumberOfCpvClusters++ ;
621 clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),0.) ; // no timing information in CPV
622 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
623 iDigitInCluster++ ;
624 fDigitsUsed[i]=kTRUE ;
625 } // else
626
627 //Now scan remaining digits in list to find neigbours of our seed
628
629 AliPHOSDigit * digitN ;
630 index = 0 ;
631 while (index < iDigitInCluster){ // scan over digits already in cluster
632 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
633 index++ ;
634 for(Int_t j=iFirst; j<nDigits; j++){
635 if (iDigitInCluster >= maxNDigits) {
636 AliError(Form("The number of digits in cluster is more than %d, skip the rest of event",
637 maxNDigits));
638 return;
639 }
640 if(fDigitsUsed[j])
641 continue ; //look through remaining digits
642 digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
643 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
644 switch (ineb ) {
645 case -1: //too early (e.g. previous module), do not look before j at subsequent passes
646 iFirst=j ;
647 break ;
648 case 0 : // not a neighbour
649 break ;
650 case 1 : // are neighbours
651 clu->AddDigit(*digitN, Calibrate(digitN->GetEnergy(),digitN->GetId()),CalibrateT(digitN->GetTime(),digitN->GetId())) ;
652 clusterdigitslist[iDigitInCluster] = j ;
653 iDigitInCluster++ ;
654 fDigitsUsed[j]=kTRUE ;
655 break ;
656 case 2 : // too far from each other
657 goto endOfLoop;
658 } // switch
659
660 }
661
662 endOfLoop: ; //scanned all possible neighbours for this digit
663
664 } // loop over cluster
665 } // energy theshold
666 }
667
668}
669
670//____________________________________________________________________________
671void AliPHOSClusterizerv1::MakeUnfolding()
672{
673 // Unfolds clusters using the shape of an ElectroMagnetic shower
674 // Performs unfolding of all EMC/CPV clusters
675
676 // Unfold first EMC clusters
677 if(fNumberOfEmcClusters > 0){
678
679 Int_t nModulesToUnfold = fGeom->GetNModules() ;
680
681 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
682 Int_t index ;
683 for(index = 0 ; index < numberofNotUnfolded ; index++){
684
685 AliPHOSEmcRecPoint * emcRecPoint = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
686 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
687 break ;
688
689 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
690 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
691 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
692 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
693
694 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
695 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
696
697 fEMCRecPoints->Remove(emcRecPoint);
698 fEMCRecPoints->Compress() ;
699 index-- ;
700 fNumberOfEmcClusters -- ;
701 numberofNotUnfolded-- ;
702 }
703 else{
704 emcRecPoint->SetNExMax(1) ; //Only one local maximum
705 }
706
707 delete[] maxAt ;
708 delete[] maxAtEnergy ;
709 }
710 }
711 // Unfolding of EMC clusters finished
712
713
714 // Unfold now CPV clusters
715 if(fNumberOfCpvClusters > 0){
716
717 Int_t nModulesToUnfold = fGeom->GetNModules() ;
718
719 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
720 Int_t index ;
721 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
722
723 AliPHOSRecPoint * recPoint = static_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
724
725 if(recPoint->GetPHOSMod()> nModulesToUnfold)
726 break ;
727
728 AliPHOSEmcRecPoint * emcRecPoint = static_cast<AliPHOSEmcRecPoint*>(recPoint) ;
729
730 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
731 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
732 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
733 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
734
735 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
736 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
737 fCPVRecPoints->Remove(emcRecPoint);
738 fCPVRecPoints->Compress() ;
739 index-- ;
740 numberofCpvNotUnfolded-- ;
741 fNumberOfCpvClusters-- ;
742 }
743
744 delete[] maxAt ;
745 delete[] maxAtEnergy ;
746 }
747 }
748 //Unfolding of Cpv clusters finished
749
750}
751
752//____________________________________________________________________________
753Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
754{
755 // Shape of the shower (see PHOS TDR)
756 // If you change this function, change also the gradient evaluation in ChiSquare()
757
758 //for the moment we neglect dependence on the incident angle.
759
760 Double_t r2 = x*x + z*z ;
761 Double_t r4 = r2*r2 ;
762 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
763 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
764 return shape ;
765}
766
767//____________________________________________________________________________
768void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
769 Int_t nMax,
770 AliPHOSDigit ** maxAt,
771 Float_t * maxAtEnergy)
772{
773 // Performs the unfolding of a cluster with nMax overlapping showers
774
775 Int_t nPar = 3 * nMax ;
776 Float_t * fitparameters = new Float_t[nPar] ;
777
778 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
779
780 if( !rv ) {
781 // Fit failed, return and remove cluster
782 iniEmc->SetNExMax(-1) ;
783 delete[] fitparameters ;
784 return ;
785 }
786
787 // create ufolded rec points and fill them with new energy lists
788 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
789 // and later correct this number in acordance with actual energy deposition
790
791 Int_t nDigits = iniEmc->GetMultiplicity() ;
792 Float_t * efit = new Float_t[nDigits] ;
793 Float_t xDigit=0.,zDigit=0. ;
794 Float_t xpar=0.,zpar=0.,epar=0. ;
795 Int_t relid[4] ;
796 AliPHOSDigit * digit = 0 ;
797 Int_t * emcDigits = iniEmc->GetDigitsList() ;
798
799 TVector3 vIncid ;
800
801 Int_t iparam ;
802 Int_t iDigit ;
803 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
804 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
805 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
806 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
807 efit[iDigit] = 0;
808
809 iparam = 0 ;
810 while(iparam < nPar ){
811 xpar = fitparameters[iparam] ;
812 zpar = fitparameters[iparam+1] ;
813 epar = fitparameters[iparam+2] ;
814 iparam += 3 ;
815// fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
816// efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
817 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
818 }
819 }
820
821 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
822 // so that energy deposited in each cell is distributed betwin new clusters proportionally
823 // to its contribution to efit
824
825 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
826 Float_t ratio ;
827
828 iparam = 0 ;
829 while(iparam < nPar ){
830 xpar = fitparameters[iparam] ;
831 zpar = fitparameters[iparam+1] ;
832 epar = fitparameters[iparam+2] ;
833 iparam += 3 ;
834// fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
835
836 AliPHOSEmcRecPoint * emcRP = 0 ;
837
838 if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
839
840 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
841 fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
842
843 (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
844 emcRP = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
845 fNumberOfEmcClusters++ ;
846 emcRP->SetNExMax((Int_t)nPar/3) ;
847 }
848 else{//create new entries in fCPVRecPoints
849 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
850 fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
851
852 (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
853 emcRP = static_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
854 fNumberOfCpvClusters++ ;
855 }
856
857 Float_t eDigit ;
858 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
859 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
860 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
861 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
862// ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
863 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
864 eDigit = emcEnergies[iDigit] * ratio ;
865 emcRP->AddDigit( *digit, eDigit,CalibrateT(digit->GetTime(),digit->GetId()) ) ;
866 }
867 }
868
869 delete[] fitparameters ;
870 delete[] efit ;
871
872}
873
874//_____________________________________________________________________________
875void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
876{
877 // Calculates the Chi square for the cluster unfolding minimization
878 // Number of parameters, Gradient, Chi squared, parameters, what to do
879
880 TList * toMinuit = static_cast<TList*>( gMinuit->GetObjectFit() ) ;
881
882 AliPHOSEmcRecPoint * emcRP = static_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
883 TClonesArray * digits = static_cast<TClonesArray*>( toMinuit->At(1) ) ;
884 // A bit buggy way to get an access to the geometry
885 // To be revised!
886 AliPHOSGeometry *geom = static_cast<AliPHOSGeometry *>(toMinuit->At(2));
887
888// TVector3 * vtx = static_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
889
890 // AliPHOSEmcRecPoint * emcRP = static_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
891
892 Int_t * emcDigits = emcRP->GetDigitsList() ;
893
894 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
895
896 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
897
898// TVector3 vInc ;
899 fret = 0. ;
900 Int_t iparam ;
901
902 if(iflag == 2)
903 for(iparam = 0 ; iparam < nPar ; iparam++)
904 Grad[iparam] = 0 ; // Will evaluate gradient
905
906 Double_t efit ;
907
908 AliPHOSDigit * digit ;
909 Int_t iDigit ;
910
911 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
912
913 digit = static_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
914
915 Int_t relid[4] ;
916 Float_t xDigit ;
917 Float_t zDigit ;
918
919 geom->AbsToRelNumbering(digit->GetId(), relid) ;
920
921 geom->RelPosInModule(relid, xDigit, zDigit) ;
922
923 if(iflag == 2){ // calculate gradient
924 Int_t iParam = 0 ;
925 efit = 0 ;
926 while(iParam < nPar ){
927 Double_t dx = (xDigit - x[iParam]) ;
928 iParam++ ;
929 Double_t dz = (zDigit - x[iParam]) ;
930 iParam++ ;
931// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
932// efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
933 efit += x[iParam] * ShowerShape(dx,dz) ;
934 iParam++ ;
935 }
936 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
937 iParam = 0 ;
938 while(iParam < nPar ){
939 Double_t xpar = x[iParam] ;
940 Double_t zpar = x[iParam+1] ;
941 Double_t epar = x[iParam+2] ;
942// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
943 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
944// Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
945 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
946//DP: No incident angle dependence in gradient yet!!!!!!
947 Double_t r4 = dr*dr*dr*dr ;
948 Double_t r295 = TMath::Power(dr,2.95) ;
949 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
950 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
951
952 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
953 iParam++ ;
954 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
955 iParam++ ;
956 Grad[iParam] += shape ; // Derivative over energy
957 iParam++ ;
958 }
959 }
960 efit = 0;
961 iparam = 0 ;
962
963 while(iparam < nPar ){
964 Double_t xpar = x[iparam] ;
965 Double_t zpar = x[iparam+1] ;
966 Double_t epar = x[iparam+2] ;
967 iparam += 3 ;
968// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
969// efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
970 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
971 }
972
973 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
974 // Here we assume, that sigma = sqrt(E)
975 }
976
977}
978
979//____________________________________________________________________________
980void AliPHOSClusterizerv1::Print(const Option_t *)const
981{
982 // Print clusterizer parameters
983
984 TString message ;
985 TString taskName(GetName()) ;
986 taskName.ReplaceAll(Version(), "") ;
987
988 if( strcmp(GetName(), "") !=0 ) {
989 // Print parameters
990 message = "\n--------------- %s %s -----------\n" ;
991 message += "Clusterizing digits from the file: %s\n" ;
992 message += " Branch: %s\n" ;
993 message += " EMC Clustering threshold = %f\n" ;
994 message += " EMC Local Maximum cut = %f\n" ;
995 message += " EMC Logarothmic weight = %f\n" ;
996 message += " CPV Clustering threshold = %f\n" ;
997 message += " CPV Local Maximum cut = %f\n" ;
998 message += " CPV Logarothmic weight = %f\n" ;
999 if(fToUnfold)
1000 message += " Unfolding on\n" ;
1001 else
1002 message += " Unfolding off\n" ;
1003
1004 message += "------------------------------------------------------------------" ;
1005 }
1006 else
1007 message = " AliPHOSClusterizerv1 not initialized " ;
1008
1009 AliInfo(Form("%s, %s %s %s %s %f %f %f %f %f %f", message.Data(),
1010 taskName.Data(),
1011 GetTitle(),
1012 taskName.Data(),
1013 GetName(),
1014 fEmcClusteringThreshold,
1015 fEmcLocMaxCut,
1016 fW0,
1017 fCpvClusteringThreshold,
1018 fCpvLocMaxCut,
1019 fW0CPV )) ;
1020}
1021//____________________________________________________________________________
1022void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1023{
1024 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1025
1026 AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
1027 fEMCRecPoints->GetEntriesFast(),
1028 fCPVRecPoints->GetEntriesFast() )) ;
1029
1030 if(strstr(option,"all")) {
1031 printf("\n EMC clusters \n") ;
1032 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1033 Int_t index ;
1034 for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1035 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1036 TVector3 locpos;
1037 rp->GetLocalPosition(locpos);
1038 Float_t lambda[2];
1039 rp->GetElipsAxis(lambda);
1040 Int_t * primaries;
1041 Int_t nprimaries;
1042 primaries = rp->GetPrimaries(nprimaries);
1043 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1044 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1045 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1046
1047 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1048 printf("%d ", primaries[iprimary] ) ;
1049 }
1050 printf("\n") ;
1051 }
1052
1053 //Now plot CPV recPoints
1054 printf("\n CPV clusters \n") ;
1055 printf("Index Ene(MeV) Module X Y Z \n") ;
1056 for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1057 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1058
1059 TVector3 locpos;
1060 rp->GetLocalPosition(locpos);
1061
1062 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1063 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1064 locpos.X(), locpos.Y(), locpos.Z()) ;
1065 }
1066 }
1067}
1068
1069
1070//____________________________________________________________________________
1071void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1072{
1073 //For each EMC rec. point set the distance to the nearest bad crystal.
1074 //Author: Boris Polichtchouk
1075
1076 if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1077
1078 Int_t badIds[8000];
1079 memset(badIds,0,8000*sizeof(Int_t));
1080 fgCalibData->EmcBadChannelIds(badIds);
1081
1082 AliPHOSEmcRecPoint* rp;
1083
1084 TMatrixF gmat;
1085 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1086 TVector3 gposBadChannel; // global position of bad crystal
1087 TVector3 dR;
1088
1089 Float_t dist,minDist;
1090 Int_t relid[4]={0,0,0,0} ;
1091 TVector3 lpos ;
1092 for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1093 rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1094 //evaluate distance to border
1095 relid[0]=rp->GetPHOSMod() ;
1096 relid[2]=1 ;
1097 relid[3]=1 ;
1098 Float_t xcorner,zcorner;
1099 fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1100 rp->GetLocalPosition(lpos) ;
1101 minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
1102 for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1103 fGeom->AbsToRelNumbering(badIds[iBad],relid) ;
1104 if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1105 continue ; //bad channels can be in the module which does not exist in simulations.
1106 rp->GetGlobalPosition(gposRecPoint,gmat);
1107 fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1108 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1109 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1110 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1111 dR = gposBadChannel-gposRecPoint;
1112 dist = dR.Mag();
1113 if(dist<minDist) minDist = dist;
1114 }
1115
1116 rp->SetDistanceToBadCrystal(minDist);
1117 }
1118
1119}
1120//==================================================================================
1121Float_t AliPHOSClusterizerv1::Calibrate(Float_t amp, Int_t absId) const{
1122 // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB
1123
1124 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1125
1126 //Determine rel.position of the cell absolute ID
1127 Int_t relId[4];
1128 geom->AbsToRelNumbering(absId,relId);
1129 Int_t module=relId[0];
1130 Int_t row =relId[2];
1131 Int_t column=relId[3];
1132 if(relId[1]){ //CPV
1133 Float_t calibration = fgCalibData->GetADCchannelCpv(module,column,row);
1134 return amp*calibration ;
1135 }
1136 else{ //EMC
1137 Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row);
1138 return amp*calibration ;
1139 }
1140}
1141//==================================================================================
1142Float_t AliPHOSClusterizerv1::CalibrateT(Float_t time, Int_t absId)const{
1143 // Calibrate time in EMC digit
1144
1145 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1146
1147 //Determine rel.position of the cell absolute ID
1148 Int_t relId[4];
1149 geom->AbsToRelNumbering(absId,relId);
1150 Int_t module=relId[0];
1151 Int_t row =relId[2];
1152 Int_t column=relId[3];
1153 if(relId[1]){ //CPV
1154 return 0. ;
1155 }
1156 else{ //EMC
1157 time += fgCalibData->GetTimeShiftEmc(module,column,row);
1158 return time ;
1159 }
1160}
1161