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