ITS DAs
[u/mrichter/AliRoot.git] / PHOS / PHOSrec / AliPHOSClusterizerv1.cxx
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
190 ClassImp(AliPHOSClusterizerv1)
191   
192 //____________________________________________________________________________
193 AliPHOSClusterizerv1::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 //____________________________________________________________________________
214 AliPHOSClusterizerv1::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 //____________________________________________________________________________
242 void 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 //____________________________________________________________________________
281 Bool_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 //____________________________________________________________________________
373 void 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 //____________________________________________________________________________
393 void 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 //____________________________________________________________________________
425 Int_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 //____________________________________________________________________________
470 Bool_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 //____________________________________________________________________________
482 Bool_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 //____________________________________________________________________________
496 Bool_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 //____________________________________________________________________________
510 void 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,fDigitsArr) ;
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 //____________________________________________________________________________
567 void 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(),digit->IsLG())) ;
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(),digit->IsLG())) ;
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 //____________________________________________________________________________
671 void 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 //____________________________________________________________________________
753 Double_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 //____________________________________________________________________________
768 void  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(),digit->IsLG()) ) ;
866     }        
867   }
868  
869   delete[] fitparameters ; 
870   delete[] efit ; 
871   
872 }
873
874 //_____________________________________________________________________________
875 void 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 //____________________________________________________________________________
980 void 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 //____________________________________________________________________________
1022 void 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 //____________________________________________________________________________
1071 void 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 //==================================================================================
1121 Float_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 //==================================================================================
1142 Float_t AliPHOSClusterizerv1::CalibrateT(Float_t time, Int_t absId,Bool_t isLG)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     if(isLG)
1158       time += fgCalibData->GetLGTimeShiftEmc(module,column,row);
1159     else
1160       time += fgCalibData->GetTimeShiftEmc(module,column,row);
1161     return time ;
1162   }
1163 }
1164