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