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