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