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