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