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