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