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