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