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