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