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