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: AliPHOSClusterizerv1.cxx,v $
  21  * Revision 1.118  2007/12/11 21:23:26  kharlov
  22  * Added possibility to swith off unfolding
  23  *
  24  * Revision 1.117  2007/10/18 08:42:05  kharlov
  25  * Bad channels cleaned before clusterization
  26  *
  27  * Revision 1.116  2007/10/01 20:24:08  kharlov
  28  * Memory leaks fixed
  29  *
  30  * Revision 1.115  2007/09/26 14:22:17  cvetan
  31  * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
  32  *
  33  * Revision 1.114  2007/09/06 16:06:44  kharlov
  34  * Absence of sorting results in loose of all unfolded clusters
  35  *
  36  * Revision 1.113  2007/08/28 12:55:07  policheh
  37  * Loaders removed from the reconstruction code (C.Cheshkov)
  38  *
  39  * Revision 1.112  2007/08/22 09:20:50  hristov
  40  * Updated QA classes (Yves)
  41  *
  42  * Revision 1.111  2007/08/08 12:11:28  kharlov
  43  * Protection against uninitialized fQADM
  44  *
  45  * Revision 1.110  2007/08/07 14:16:00  kharlov
  46  * Quality assurance added (Yves Schutz)
  47  *
  48  * Revision 1.109  2007/07/24 17:20:35  policheh
  49  * Usage of RecoParam objects instead of hardcoded parameters in reconstruction.
  50  * (See $ALICE_ROOT/PHOS/macros/BeamTest2006/RawReconstruction.C).
  51  *
  52  * Revision 1.108  2007/06/18 07:00:51  kharlov
  53  * Bug fix for attempt to use AliPHOSEmcRecPoint after its deletion
  54  *
  55  * Revision 1.107  2007/05/25 14:12:26  policheh
  56  * Local to tracking CS transformation added for CPV rec. points
  57  *
  58  * Revision 1.106  2007/05/24 13:01:22  policheh
  59  * Local to tracking CS transformation invoked for each EMC rec.point
  60  *
  61  * Revision 1.105  2007/05/02 13:41:22  kharlov
  62  * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
  63  *
  64  * Revision 1.104  2007/04/27 16:55:53  kharlov
  65  * Calibration stops if PHOS CDB objects do not exist
  66  *
  67  * Revision 1.103  2007/04/11 11:55:45  policheh
  68  * SetDistancesToBadChannels() added.
  69  *
  70  * Revision 1.102  2007/03/28 19:18:15  kharlov
  71  * RecPoints recalculation in TSM removed
  72  *
  73  * Revision 1.101  2007/03/06 06:51:27  kharlov
  74  * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
  75  *
  76  * Revision 1.100  2007/01/10 11:07:26  kharlov
  77  * Raw digits writing to file (B.Polichtchouk)
  78  *
  79  * Revision 1.99  2006/11/07 16:49:51  kharlov
  80  * Corrections for next event switch in case of raw data (B.Polichtchouk)
  81  *
  82  * Revision 1.98  2006/10/27 17:14:27  kharlov
  83  * Introduce AliDebug and AliLog (B.Polichtchouk)
  84  *
  85  * Revision 1.97  2006/08/29 11:41:19  kharlov
  86  * Missing implementation of ctors and = operator are added
  87  *
  88  * Revision 1.96  2006/08/25 16:56:30  kharlov
  89  * Compliance with Effective C++
  90  *
  91  * Revision 1.95  2006/08/11 12:36:26  cvetan
  92  * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
  93  *
  94  * Revision 1.94  2006/08/07 12:27:49  hristov
  95  * Removing obsolete code which affected the event numbering scheme
  96  *
  97  * Revision 1.93  2006/08/01 12:20:17  cvetan
  98  * 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
  99  *
 100  * Revision 1.92  2006/04/29 20:26:46  hristov
 101  * Separate EMC and CPV calibration (Yu.Kharlov)
 102  *
 103  * Revision 1.91  2006/04/22 10:30:17  hristov
 104  * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
 105  *
 106  * Revision 1.90  2006/04/11 15:22:59  hristov
 107  * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
 108  *
 109  * Revision 1.89  2006/03/13 14:05:42  kharlov
 110  * Calibration objects for EMC and CPV
 111  *
 112  * Revision 1.88  2006/01/11 08:54:52  hristov
 113  * Additional protection in case no calibration entry was found
 114  *
 115  * Revision 1.87  2005/11/22 08:46:43  kharlov
 116  * Updated to new CDB (Boris Polichtchouk)
 117  *
 118  * Revision 1.86  2005/11/14 21:52:43  hristov
 119  * Coding conventions
 120  *
 121  * Revision 1.85  2005/09/27 16:08:08  hristov
 122  * New version of CDB storage framework (A.Colla)
 123  *
 124  * Revision 1.84  2005/09/21 10:02:47  kharlov
 125  * Reading calibration from CDB (Boris Polichtchouk)
 126  *
 127  * Revision 1.82  2005/09/02 15:43:13  kharlov
 128  * Add comments in GetCalibrationParameters and Calibrate
 129  *
 130  * Revision 1.81  2005/09/02 14:32:07  kharlov
 131  * Calibration of raw data
 132  *
 133  * Revision 1.80  2005/08/24 15:31:36  kharlov
 134  * Setting raw digits flag
 135  *
 136  * Revision 1.79  2005/07/25 15:53:53  kharlov
 137  * Read raw data
 138  *
 139  * Revision 1.78  2005/05/28 14:19:04  schutz
 140  * Compilation warnings fixed by T.P.
 141  *
 142  */
 143
 144 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
 145 //////////////////////////////////////////////////////////////////////////////
 146 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and
 147 //  unfolds the clusters having several local maxima.
 148 //  Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
 149 //  PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
 150 //  parameters including input digits branch title, thresholds etc.)
 151 //  This TTask is normally called from Reconstructioner, but can as well be used in
 152 //  standalone mode.
 153 // Use Case:
 154 //  root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1(<pointer_to_phos_geometry_onject>)
 155 //  root [1] cl->Digits2Clusters(digitsTree,clusterTree)
 156 //               //finds RecPoints in the current event
 157 //  root [2] cl->SetDigitsBranch("digits2")
 158 //               //sets another title for Digitis (input) branch
 159 //  root [3] cl->SetRecPointsBranch("recp2")
 160 //               //sets another title four output branches
 161 //  root [4] cl->SetEmcLocalMaxCut(0.03)
 162 //               //set clusterization parameters
 163
 164 // --- ROOT system ---
 165
 166 #include "TMath.h"
 167 #include "TMinuit.h"
 168 #include "TTree.h"
 169 #include "TBenchmark.h"
 170 #include "TClonesArray.h"
 171
 172 // --- Standard library ---
 173
 174 // --- AliRoot header files ---
 175 #include "AliLog.h"
 176 #include "AliConfig.h"
 177 #include "AliPHOSGeometry.h"
 178 #include "AliPHOSClusterizerv1.h"
 179 #include "AliPHOSEmcRecPoint.h"
 180 #include "AliPHOSCpvRecPoint.h"
 181 #include "AliPHOSDigit.h"
 182 #include "AliPHOSDigitizer.h"
 183 #include "AliPHOSCalibrationDB.h"
 184 #include "AliCDBManager.h"
 185 #include "AliCDBStorage.h"
 186 #include "AliCDBEntry.h"
 187 #include "AliPHOSRecoParam.h"
 188 #include "AliPHOSReconstructor.h"
 189 #include "AliPHOSCalibData.h"
 190
 191 ClassImp(AliPHOSClusterizerv1)
 192
 193 //____________________________________________________________________________
 194 AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
 195   AliPHOSClusterizer(),
 196   fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
 197   fWrite(0),
 198   fNumberOfEmcClusters(0),    fNumberOfCpvClusters(0),
 199   fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
 200   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
 201   fW0CPV(0),                  fEmcTimeGate(0),          fEcoreRadius(0)
 202 {
 203   // default ctor (to be used mainly by Streamer)
 204
 205   fDefaultInit = kTRUE ;
 206 }
 207
 208 //____________________________________________________________________________
 209 AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
 210   AliPHOSClusterizer(geom),
 211   fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
 212   fWrite(0),
 213   fNumberOfEmcClusters(0),    fNumberOfCpvClusters(0),
 214   fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
 215   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
 216   fW0CPV(0),                  fEmcTimeGate(0),          fEcoreRadius(0)
 217 {
 218   // ctor with the indication of the file where header Tree and digits Tree are stored
 219
 220   Init() ;
 221   fDefaultInit = kFALSE ;
 222 }
 223
 224 //____________________________________________________________________________
 225   AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
 226 {
 227   // dtor
 228
 229 }
 230 //____________________________________________________________________________
 231 void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
 232 {
 233   // Steering method to perform clusterization for one event
 234   // The input is the tree with digits.
 235   // The output is the tree with clusters.
 236
 237   if(strstr(option,"tim"))
 238     gBenchmark->Start("PHOSClusterizer");
 239
 240   if(strstr(option,"print")) {
 241     Print() ;
 242     return ;
 243   }
 244
 245   MakeClusters() ;
 246
 247   AliDebug(2,Form(" ---- Printing clusters (%d)\n",
 248                   fEMCRecPoints->GetEntries()));
 249   if(AliLog::GetGlobalDebugLevel()>1)
 250     fEMCRecPoints->Print();
 251
 252   if(fToUnfold)
 253     MakeUnfolding();
 254
 255   WriteRecPoints();
 256
 257   if(strstr(option,"deb"))
 258     PrintRecPoints(option) ;
 259
 260   if(strstr(option,"tim")){
 261     gBenchmark->Stop("PHOSClusterizer");
 262     AliInfo(Form("took %f seconds for Clusterizing\n",
 263                  gBenchmark->GetCpuTime("PHOSClusterizer")));
 264   }
 265   fEMCRecPoints->Delete();
 266   fCPVRecPoints->Delete();
 267 }
 268
 269 //____________________________________________________________________________
 270 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
 271                                     Int_t nPar, Float_t * fitparameters) const
 272 {
 273   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 274   // The initial values for fitting procedure are set equal to the positions of local maxima.
 275   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 276
 277
 278   gMinuit->mncler();                     // Reset Minuit's list of paramters
 279   gMinuit->SetPrintLevel(-1) ;           // No Printout
 280   gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
 281                                          // To set the address of the minimization function
 282
 283   TList * toMinuit = new TList();
 284   toMinuit->AddAt(emcRP,0) ;
 285   toMinuit->AddAt(fDigitsArr,1) ;
 286   toMinuit->AddAt(fGeom,2) ;
 287
 288   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 289
 290   // filling initial values for fit parameters
 291   AliPHOSDigit * digit ;
 292
 293   Int_t ierflg  = 0;
 294   Int_t index   = 0 ;
 295   Int_t nDigits = (Int_t) nPar / 3 ;
 296
 297   Int_t iDigit ;
 298
 299   for(iDigit = 0; iDigit < nDigits; iDigit++){
 300     digit = maxAt[iDigit];
 301
 302     Int_t relid[4] ;
 303     Float_t x = 0.;
 304     Float_t z = 0.;
 305     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 306     fGeom->RelPosInModule(relid, x, z) ;
 307
 308     Float_t energy = maxAtEnergy[iDigit] ;
 309
 310     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 311     index++ ;
 312     if(ierflg != 0){
 313       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
 314       return kFALSE;
 315     }
 316     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 317     index++ ;
 318     if(ierflg != 0){
 319        Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
 320       return kFALSE;
 321     }
 322     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 323     index++ ;
 324     if(ierflg != 0){
 325       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
 326       return kFALSE;
 327     }
 328   }
 329
 330   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 331                       //  depends on it.
 332   Double_t p1 = 1.0 ;
 333   Double_t p2 = 0.0 ;
 334
 335   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 336   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 337   gMinuit->SetMaxIterations(5);
 338   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 339
 340   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 341
 342   if(ierflg == 4){  // Minimum not found
 343     Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
 344     return kFALSE ;
 345   }
 346   for(index = 0; index < nPar; index++){
 347     Double_t err ;
 348     Double_t val ;
 349     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 350     fitparameters[index] = val ;
 351    }
 352
 353   delete toMinuit ;
 354   return kTRUE;
 355
 356 }
 357
 358
 359 //____________________________________________________________________________
 360 void AliPHOSClusterizerv1::Init()
 361 {
 362   // Make all memory allocations which can not be done in default constructor.
 363   // Attach the Clusterizer task to the list of PHOS tasks
 364
 365   fEmcCrystals = fGeom->GetNModules() *  fGeom->GetNCristalsInModule() ;
 366
 367   if(!gMinuit)
 368     gMinuit = new TMinuit(100);
 369
 370   if (!fgCalibData)
 371     fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
 372   if (fgCalibData->GetCalibDataEmc() == 0)
 373     AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
 374
 375 }
 376
 377 //____________________________________________________________________________
 378 void AliPHOSClusterizerv1::InitParameters()
 379 {
 380
 381   fNumberOfCpvClusters     = 0 ;
 382   fNumberOfEmcClusters     = 0 ;
 383
 384   const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
 385   if(!recoParam) AliFatal("Reconstruction parameters are not set!");
 386
 387   recoParam->Print();
 388
 389   fEmcClusteringThreshold  = recoParam->GetEMCClusteringThreshold();
 390   fCpvClusteringThreshold  = recoParam->GetCPVClusteringThreshold();
 391
 392   fEmcLocMaxCut            = recoParam->GetEMCLocalMaxCut();
 393   fCpvLocMaxCut            = recoParam->GetCPVLocalMaxCut();
 394
 395   fW0                      = recoParam->GetEMCLogWeight();
 396   fW0CPV                   = recoParam->GetCPVLogWeight();
 397
 398   fEmcTimeGate             = 1.e-6 ;
 399   fEcoreRadius             = recoParam->GetEMCEcoreRadius();
 400
 401   fToUnfold                = recoParam->EMCToUnfold() ;
 402
 403   fWrite                   = kTRUE ;
 404 }
 405
 406 //____________________________________________________________________________
 407 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
 408 {
 409   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 410   //                                       = 1 are neighbour
 411   //                                       = 2 are not neighbour but do not continue searching
 412   //                                       =-1 are not neighbour, continue searching, but do not look before d2 next time
 413   // neighbours are defined as digits having at least a common vertex
 414   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 415   //                                      which is compared to a digit (d2)  not yet in a cluster
 416
 417   Int_t relid1[4] ;
 418   fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
 419
 420   Int_t relid2[4] ;
 421   fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
 422
 423   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
 424     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 425     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 426
 427     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){   //At least common vertex
 428       //    if (( relid1[2]==relid2[2] && coldiff <= 1 )  || ( relid1[3]==relid2[3] &&  rowdiff <= 1 )){ //common side
 429       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
 430       return 1 ;
 431     }
 432     else {
 433       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
 434         return 2; //  Difference in row numbers is too large to look further
 435     }
 436     return 0 ;
 437
 438   }
 439   else {
 440     if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
 441       return -1 ;
 442     if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
 443       return -1 ;
 444
 445     return 2 ;
 446
 447   }
 448
 449   return 0 ;
 450 }
 451 //____________________________________________________________________________
 452 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
 453 {
 454   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
 455
 456   Bool_t rv = kFALSE ;
 457
 458   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 459
 460   if(digit->GetId() <= nEMC )   rv = kTRUE;
 461
 462   return rv ;
 463 }
 464
 465 //____________________________________________________________________________
 466 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 467 {
 468   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
 469
 470   Bool_t rv = kFALSE ;
 471
 472   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 473
 474   if(digit->GetId() > nEMC )   rv = kTRUE;
 475
 476   return rv ;
 477 }
 478
 479 //____________________________________________________________________________
 480 void AliPHOSClusterizerv1::WriteRecPoints()
 481 {
 482
 483   // Creates new branches with given title
 484   // fills and writes into TreeR.
 485
 486   Int_t index ;
 487   //Evaluate position, dispersion and other RecPoint properties..
 488   Int_t nEmc = fEMCRecPoints->GetEntriesFast();
 489   Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
 490   TVector3 fakeVtx(0.,0.,0.) ;
 491   for(index = 0; index < nEmc; index++){
 492     AliPHOSEmcRecPoint * rp =
 493       dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
 494     rp->Purify(emcMinE) ;
 495     if(rp->GetMultiplicity()==0){
 496       fEMCRecPoints->RemoveAt(index) ;
 497       delete rp ;
 498       continue;
 499     }
 500
 501     // No vertex is available now, calculate corrections in PID
 502     rp->EvalAll(fDigitsArr) ;
 503     rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
 504     rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
 505     rp->EvalLocal2TrackingCSTransform();
 506   }
 507   fEMCRecPoints->Compress() ;
 508   fEMCRecPoints->Sort() ;
 509   //  fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
 510   for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
 511     dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
 512   }
 513
 514   //For each rec.point set the distance to the nearest bad crystal (BVP)
 515   SetDistancesToBadChannels();
 516
 517   //Now the same for CPV
 518   for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
 519     AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
 520     rp->EvalAll(fDigitsArr) ;
 521     rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
 522     rp->EvalLocal2TrackingCSTransform();
 523   }
 524   fCPVRecPoints->Sort() ;
 525
 526   for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
 527     dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
 528
 529   fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
 530
 531   if(fWrite){ //We write TreeR
 532     fTreeR->Fill();
 533   }
 534 }
 535
 536 //____________________________________________________________________________
 537 void AliPHOSClusterizerv1::MakeClusters()
 538 {
 539   // Steering method to construct the clusters stored in a list of Reconstructed Points
 540   // A cluster is defined as a list of neighbour digits
 541
 542   fNumberOfCpvClusters     = 0 ;
 543   fNumberOfEmcClusters     = 0 ;
 544
 545   //Mark all digits as unused yet
 546   const Int_t maxNDigits = 1500;
 547   Int_t nDigits=fDigitsArr->GetEntriesFast() ;
 548   if (nDigits > maxNDigits) {
 549     AliWarning(Form("Skip event with too high digit occupancy: nDigits=%d",nDigits));
 550     return;
 551   }
 552
 553   for(Int_t i=0; i<nDigits; i++){
 554     fDigitsUsed[i]=0 ;
 555   }
 556   Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
 557                      //e.g. first digit in this module, first CPV digit etc.
 558   AliPHOSDigit * digit ;
 559   TArrayI clusterdigitslist(maxNDigits) ;
 560   AliPHOSRecPoint * clu = 0 ;
 561   for(Int_t i=0; i<nDigits; i++){
 562     if(fDigitsUsed[i])
 563       continue ;
 564
 565     digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
 566
 567     clu=0 ;
 568
 569     Int_t index ;
 570
 571     //is this digit so energetic that start cluster?
 572     if (( IsInEmc(digit) &&  digit->GetEnergy() > fEmcClusteringThreshold ) ||
 573         ( IsInCpv(digit) &&  digit->GetEnergy() > fCpvClusteringThreshold ) ) {
 574       Int_t iDigitInCluster = 0 ;
 575
 576       if  ( IsInEmc(digit) ) {
 577         // start a new EMC RecPoint
 578         if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 579           fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
 580
 581         fEMCRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
 582         clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 583         fNumberOfEmcClusters++ ;
 584         clu->AddDigit(*digit, digit->GetEnergy()) ;
 585         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
 586         iDigitInCluster++ ;
 587         fDigitsUsed[i]=kTRUE ;
 588       } else {
 589         // start a new CPV cluster
 590         if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 591           fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
 592
 593         fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
 594         clu =  static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 595         fNumberOfCpvClusters++ ;
 596         clu->AddDigit(*digit, digit->GetEnergy())  ;
 597         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;
 598         iDigitInCluster++ ;
 599         fDigitsUsed[i]=kTRUE ;
 600       } // else
 601
 602       //Now scan remaining digits in list to find neigbours of our seed
 603
 604       AliPHOSDigit * digitN ;
 605       index = 0 ;
 606       while (index < iDigitInCluster){ // scan over digits already in cluster
 607         digit =  static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) )  ;
 608         index++ ;
 609         for(Int_t j=iFirst; j<nDigits; j++){
 610           if(fDigitsUsed[j])
 611             continue ;        //look through remaining digits
 612           digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
 613           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 614           switch (ineb ) {
 615           case -1:   //too early (e.g. previous module), do not look before j at subsequent passes
 616             iFirst=j ;
 617             break ;
 618           case 0 :   // not a neighbour
 619             break ;
 620           case 1 :   // are neighbours
 621             clu->AddDigit(*digitN, digitN->GetEnergy());
 622             clusterdigitslist[iDigitInCluster] = j ;
 623             iDigitInCluster++ ;
 624             fDigitsUsed[j]=kTRUE ;
 625             break ;
 626           case 2 :   // too far from each other
 627             goto endofloop;
 628           } // switch
 629
 630         }
 631
 632         endofloop: ; //scanned all possible neighbours for this digit
 633
 634       } // loop over cluster
 635     } // energy theshold
 636   }
 637
 638 }
 639
 640 //____________________________________________________________________________
 641 void AliPHOSClusterizerv1::MakeUnfolding()
 642 {
 643   // Unfolds clusters using the shape of an ElectroMagnetic shower
 644   // Performs unfolding of all EMC/CPV clusters
 645
 646   // Unfold first EMC clusters
 647   if(fNumberOfEmcClusters > 0){
 648
 649     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 650
 651     Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
 652     Int_t index ;
 653     for(index = 0 ; index < numberofNotUnfolded ; index++){
 654
 655       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
 656       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
 657         break ;
 658
 659       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 660       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 661       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 662       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
 663
 664       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 665         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 666
 667         fEMCRecPoints->Remove(emcRecPoint);
 668         fEMCRecPoints->Compress() ;
 669         index-- ;
 670         fNumberOfEmcClusters -- ;
 671         numberofNotUnfolded-- ;
 672       }
 673       else{
 674         emcRecPoint->SetNExMax(1) ; //Only one local maximum
 675       }
 676
 677       delete[] maxAt ;
 678       delete[] maxAtEnergy ;
 679     }
 680   }
 681   // Unfolding of EMC clusters finished
 682
 683
 684   // Unfold now CPV clusters
 685   if(fNumberOfCpvClusters > 0){
 686
 687     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 688
 689     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
 690     Int_t index ;
 691     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
 692
 693       AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
 694
 695       if(recPoint->GetPHOSMod()> nModulesToUnfold)
 696         break ;
 697
 698       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
 699
 700       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 701       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 702       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 703       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
 704
 705       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 706         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 707         fCPVRecPoints->Remove(emcRecPoint);
 708         fCPVRecPoints->Compress() ;
 709         index-- ;
 710         numberofCpvNotUnfolded-- ;
 711         fNumberOfCpvClusters-- ;
 712       }
 713
 714       delete[] maxAt ;
 715       delete[] maxAtEnergy ;
 716     }
 717   }
 718   //Unfolding of Cpv clusters finished
 719
 720 }
 721
 722 //____________________________________________________________________________
 723 Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
 724 {
 725   // Shape of the shower (see PHOS TDR)
 726   // If you change this function, change also the gradient evaluation in ChiSquare()
 727
 728   //for the moment we neglect dependence on the incident angle.
 729
 730   Double_t r2    = x*x + z*z ;
 731   Double_t r4    = r2*r2 ;
 732   Double_t r295  = TMath::Power(r2, 2.95/2.) ;
 733   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 734   return shape ;
 735 }
 736
 737 //____________________________________________________________________________
 738 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 739                                           Int_t nMax,
 740                                           AliPHOSDigit ** maxAt,
 741                                           Float_t * maxAtEnergy)
 742 {
 743   // Performs the unfolding of a cluster with nMax overlapping showers
 744
 745   Int_t nPar = 3 * nMax ;
 746   Float_t * fitparameters = new Float_t[nPar] ;
 747
 748   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
 749
 750   if( !rv ) {
 751     // Fit failed, return and remove cluster
 752     iniEmc->SetNExMax(-1) ;
 753     delete[] fitparameters ;
 754     return ;
 755   }
 756
 757   // create ufolded rec points and fill them with new energy lists
 758   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
 759   // and later correct this number in acordance with actual energy deposition
 760
 761   Int_t nDigits = iniEmc->GetMultiplicity() ;
 762   Float_t * efit = new Float_t[nDigits] ;
 763   Float_t xDigit=0.,zDigit=0. ;
 764   Float_t xpar=0.,zpar=0.,epar=0.  ;
 765   Int_t relid[4] ;
 766   AliPHOSDigit * digit = 0 ;
 767   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 768
 769   TVector3 vIncid ;
 770
 771   Int_t iparam ;
 772   Int_t iDigit ;
 773   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 774     digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
 775     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 776     fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 777     efit[iDigit] = 0;
 778
 779     iparam = 0 ;
 780     while(iparam < nPar ){
 781       xpar = fitparameters[iparam] ;
 782       zpar = fitparameters[iparam+1] ;
 783       epar = fitparameters[iparam+2] ;
 784       iparam += 3 ;
 785 //      fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
 786 //      efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
 787       efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 788     }
 789   }
 790
 791   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
 792   // so that energy deposited in each cell is distributed betwin new clusters proportionally
 793   // to its contribution to efit
 794
 795   Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
 796   Float_t ratio ;
 797
 798   iparam = 0 ;
 799   while(iparam < nPar ){
 800     xpar = fitparameters[iparam] ;
 801     zpar = fitparameters[iparam+1] ;
 802     epar = fitparameters[iparam+2] ;
 803     iparam += 3 ;
 804 //    fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
 805
 806     AliPHOSEmcRecPoint * emcRP = 0 ;
 807
 808     if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
 809
 810       if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 811         fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
 812
 813       (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
 814       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 815       fNumberOfEmcClusters++ ;
 816       emcRP->SetNExMax((Int_t)nPar/3) ;
 817     }
 818     else{//create new entries in fCPVRecPoints
 819       if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 820         fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
 821
 822       (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
 823       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 824       fNumberOfCpvClusters++ ;
 825     }
 826
 827     Float_t eDigit ;
 828     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 829       digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
 830       fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 831       fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 832 //      ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
 833       ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
 834       eDigit = emcEnergies[iDigit] * ratio ;
 835       emcRP->AddDigit( *digit, eDigit ) ;
 836     }
 837   }
 838
 839   delete[] fitparameters ;
 840   delete[] efit ;
 841
 842 }
 843
 844 //_____________________________________________________________________________
 845 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 846 {
 847   // Calculates the Chi square for the cluster unfolding minimization
 848   // Number of parameters, Gradient, Chi squared, parameters, what to do
 849
 850   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
 851
 852   AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
 853   TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
 854   // A bit buggy way to get an access to the geometry
 855   // To be revised!
 856   AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
 857
 858 //  TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ;  //Vertex position
 859
 860   //  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
 861
 862   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 863
 864   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
 865
 866   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
 867
 868 //  TVector3 vInc ;
 869   fret = 0. ;
 870   Int_t iparam ;
 871
 872   if(iflag == 2)
 873     for(iparam = 0 ; iparam < nPar ; iparam++)
 874       Grad[iparam] = 0 ; // Will evaluate gradient
 875
 876   Double_t efit ;
 877
 878   AliPHOSDigit * digit ;
 879   Int_t iDigit ;
 880
 881   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 882
 883     digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
 884
 885     Int_t relid[4] ;
 886     Float_t xDigit ;
 887     Float_t zDigit ;
 888
 889     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 890
 891     geom->RelPosInModule(relid, xDigit, zDigit) ;
 892
 893      if(iflag == 2){  // calculate gradient
 894        Int_t iParam = 0 ;
 895        efit = 0 ;
 896        while(iParam < nPar ){
 897          Double_t dx = (xDigit - x[iParam]) ;
 898          iParam++ ;
 899          Double_t dz = (zDigit - x[iParam]) ;
 900          iParam++ ;
 901 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
 902 //         efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
 903          efit += x[iParam] * ShowerShape(dx,dz) ;
 904          iParam++ ;
 905        }
 906        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
 907        iParam = 0 ;
 908        while(iParam < nPar ){
 909          Double_t xpar = x[iParam] ;
 910          Double_t zpar = x[iParam+1] ;
 911          Double_t epar = x[iParam+2] ;
 912 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
 913          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
 914 //         Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
 915          Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 916 //DP: No incident angle dependence in gradient yet!!!!!!
 917          Double_t r4 = dr*dr*dr*dr ;
 918          Double_t r295 = TMath::Power(dr,2.95) ;
 919          Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
 920                                          0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
 921
 922          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
 923          iParam++ ;
 924          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
 925          iParam++ ;
 926          Grad[iParam] += shape ;                                  // Derivative over energy
 927          iParam++ ;
 928        }
 929      }
 930      efit = 0;
 931      iparam = 0 ;
 932
 933      while(iparam < nPar ){
 934        Double_t xpar = x[iparam] ;
 935        Double_t zpar = x[iparam+1] ;
 936        Double_t epar = x[iparam+2] ;
 937        iparam += 3 ;
 938 //       fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
 939 //       efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
 940        efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 941      }
 942
 943      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
 944      // Here we assume, that sigma = sqrt(E)
 945   }
 946
 947 }
 948
 949 //____________________________________________________________________________
 950 void AliPHOSClusterizerv1::Print(const Option_t *)const
 951 {
 952   // Print clusterizer parameters
 953
 954   TString message ;
 955   TString taskName(GetName()) ;
 956   taskName.ReplaceAll(Version(), "") ;
 957
 958   if( strcmp(GetName(), "") !=0 ) {
 959     // Print parameters
 960     message  = "\n--------------- %s %s -----------\n" ;
 961     message += "Clusterizing digits from the file: %s\n" ;
 962     message += "                           Branch: %s\n" ;
 963     message += "                       EMC Clustering threshold = %f\n" ;
 964     message += "                       EMC Local Maximum cut    = %f\n" ;
 965     message += "                       EMC Logarothmic weight   = %f\n" ;
 966     message += "                       CPV Clustering threshold = %f\n" ;
 967     message += "                       CPV Local Maximum cut    = %f\n" ;
 968     message += "                       CPV Logarothmic weight   = %f\n" ;
 969     if(fToUnfold)
 970       message += " Unfolding on\n" ;
 971     else
 972       message += " Unfolding off\n" ;
 973
 974     message += "------------------------------------------------------------------" ;
 975   }
 976   else
 977     message = " AliPHOSClusterizerv1 not initialized " ;
 978
 979   AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
 980        taskName.Data(),
 981        GetTitle(),
 982        taskName.Data(),
 983        GetName(),
 984        fEmcClusteringThreshold,
 985        fEmcLocMaxCut,
 986        fW0,
 987        fCpvClusteringThreshold,
 988        fCpvLocMaxCut,
 989        fW0CPV )) ;
 990 }
 991 //____________________________________________________________________________
 992 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
 993 {
 994   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
 995
 996   AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
 997                fEMCRecPoints->GetEntriesFast(),
 998                fCPVRecPoints->GetEntriesFast() ))  ;
 999
1000   if(strstr(option,"all")) {
1001     printf("\n EMC clusters \n") ;
1002     printf("Index    Ene(MeV) Multi Module     X    Y   Z    Lambdas_1  Lambda_2  # of prim  Primaries list\n") ;
1003     Int_t index ;
1004     for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1005       AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1006       TVector3  locpos;
1007       rp->GetLocalPosition(locpos);
1008       Float_t lambda[2];
1009       rp->GetElipsAxis(lambda);
1010       Int_t * primaries;
1011       Int_t nprimaries;
1012       primaries = rp->GetPrimaries(nprimaries);
1013       printf("\n%6d  %8.2f  %3d     %2d     %4.1f   %4.1f   %4.1f   %4f  %4f    %2d     : ",
1014               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1015               locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1016
1017       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1018         printf("%d ", primaries[iprimary] ) ;
1019       }
1020       printf("\n") ;
1021     }
1022
1023     //Now plot CPV recPoints
1024     printf("\n CPV clusters \n") ;
1025     printf("Index    Ene(MeV) Module     X     Y    Z  \n") ;
1026     for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1027       AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1028
1029       TVector3  locpos;
1030       rp->GetLocalPosition(locpos);
1031
1032       printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f \n",
1033              rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1034              locpos.X(), locpos.Y(), locpos.Z()) ;
1035     }
1036   }
1037 }
1038
1039
1040 //____________________________________________________________________________
1041 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1042 {
1043   //For each EMC rec. point set the distance to the nearest bad crystal.
1044   //Author: Boris Polichtchouk
1045
1046   if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1047   AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
1048
1049   Int_t badIds[8000];
1050   fgCalibData->EmcBadChannelIds(badIds);
1051
1052   AliPHOSEmcRecPoint* rp;
1053
1054   TMatrixF gmat;
1055   TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1056   TVector3 gposBadChannel; // global position of bad crystal
1057   TVector3 dR;
1058
1059   Float_t dist,minDist;
1060   Int_t relid[4]={0,0,0,0} ;
1061   TVector3 lpos ;
1062   for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1063     rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1064     //evaluate distance to border
1065     relid[0]=rp->GetPHOSMod() ;
1066     relid[2]=1 ;
1067     relid[3]=1 ;
1068     Float_t xcorner,zcorner;
1069     fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1070     rp->GetLocalPosition(lpos) ;
1071     minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
1072     for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1073       fGeom->AbsToRelNumbering(badIds[iBad],relid)  ;
1074       if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1075         continue ;                   //bad channels can be in the module which does not exist in simulations.
1076       rp->GetGlobalPosition(gposRecPoint,gmat);
1077       fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1078       AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1079                       gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1080                       gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1081       dR = gposBadChannel-gposRecPoint;
1082       dist = dR.Mag();
1083       if(dist<minDist) minDist = dist;
1084     }
1085
1086     rp->SetDistanceToBadCrystal(minDist);
1087   }
1088
1089 }