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