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