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