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