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