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