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