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   fCpvClusteringThreshold  = recoParam->GetEMCClusteringThreshold();
 388   fEmcClusteringThreshold  = recoParam->GetCPVClusteringThreshold();
 389
 390   fEmcLocMaxCut            = recoParam->GetEMCLocalMaxCut();
 391   fCpvLocMaxCut            = recoParam->GetCPVLocalMaxCut();
 392
 393   fW0                      = recoParam->GetEMCLogWeight();
 394   fW0CPV                   = recoParam->GetCPVLogWeight();
 395
 396   fEmcTimeGate             = 1.e-6 ;
 397   fEcoreRadius             = recoParam->GetEMCEcoreRadius();
 398
 399   fToUnfold                = recoParam->EMCToUnfold() ;
 400
 401   fWrite                   = kTRUE ;
 402 }
 403
 404 //____________________________________________________________________________
 405 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
 406 {
 407   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 408   //                                       = 1 are neighbour
 409   //                                       = 2 are not neighbour but do not continue searching
 410   //                                       =-1 are not neighbour, continue searching, but do not look before d2 next time
 411   // neighbours are defined as digits having at least a common vertex
 412   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 413   //                                      which is compared to a digit (d2)  not yet in a cluster
 414
 415   Int_t relid1[4] ;
 416   fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
 417
 418   Int_t relid2[4] ;
 419   fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
 420
 421   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
 422     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 423     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 424
 425     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){   //At least common vertex
 426       //    if (( relid1[2]==relid2[2] && coldiff <= 1 )  || ( relid1[3]==relid2[3] &&  rowdiff <= 1 )){ //common side
 427       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
 428       return 1 ;
 429     }
 430     else {
 431       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
 432         return 2; //  Difference in row numbers is too large to look further
 433     }
 434     return 0 ;
 435
 436   }
 437   else {
 438     if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
 439       return -1 ;
 440     if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
 441       return -1 ;
 442
 443     return 2 ;
 444
 445   }
 446
 447   return 0 ;
 448 }
 449 //____________________________________________________________________________
 450 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
 451 {
 452   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
 453
 454   Bool_t rv = kFALSE ;
 455
 456   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 457
 458   if(digit->GetId() <= nEMC )   rv = kTRUE;
 459
 460   return rv ;
 461 }
 462
 463 //____________________________________________________________________________
 464 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 465 {
 466   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
 467
 468   Bool_t rv = kFALSE ;
 469
 470   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 471
 472   if(digit->GetId() > nEMC )   rv = kTRUE;
 473
 474   return rv ;
 475 }
 476
 477 //____________________________________________________________________________
 478 void AliPHOSClusterizerv1::WriteRecPoints()
 479 {
 480
 481   // Creates new branches with given title
 482   // fills and writes into TreeR.
 483
 484   Int_t index ;
 485   //Evaluate position, dispersion and other RecPoint properties..
 486   Int_t nEmc = fEMCRecPoints->GetEntriesFast();
 487   Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
 488   TVector3 fakeVtx(0.,0.,0.) ;
 489   for(index = 0; index < nEmc; index++){
 490     AliPHOSEmcRecPoint * rp =
 491       dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
 492     rp->Purify(emcMinE) ;
 493     if(rp->GetMultiplicity()==0){
 494       fEMCRecPoints->RemoveAt(index) ;
 495       delete rp ;
 496       continue;
 497     }
 498
 499     // No vertex is available now, calculate corrections in PID
 500     rp->EvalAll(fDigitsArr) ;
 501     rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
 502     rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
 503     rp->EvalLocal2TrackingCSTransform();
 504   }
 505   fEMCRecPoints->Compress() ;
 506   fEMCRecPoints->Sort() ;
 507   //  fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
 508   for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
 509     dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
 510   }
 511
 512   //For each rec.point set the distance to the nearest bad crystal (BVP)
 513   SetDistancesToBadChannels();
 514
 515   //Now the same for CPV
 516   for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
 517     AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
 518     rp->EvalAll(fDigitsArr) ;
 519     rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
 520     rp->EvalLocal2TrackingCSTransform();
 521   }
 522   fCPVRecPoints->Sort() ;
 523
 524   for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
 525     dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
 526
 527   fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
 528
 529   if(fWrite){ //We write TreeR
 530     fTreeR->Fill();
 531   }
 532 }
 533
 534 //____________________________________________________________________________
 535 void AliPHOSClusterizerv1::MakeClusters()
 536 {
 537   // Steering method to construct the clusters stored in a list of Reconstructed Points
 538   // A cluster is defined as a list of neighbour digits
 539
 540   fNumberOfCpvClusters     = 0 ;
 541   fNumberOfEmcClusters     = 0 ;
 542
 543   //Mark all digits as unused yet
 544   const Int_t maxNDigits = 1500;
 545   Int_t nDigits=fDigitsArr->GetEntriesFast() ;
 546   if (nDigits > maxNDigits) {
 547     AliWarning(Form("Skip event with too high digit occupancy: nDigits=%d",nDigits));
 548     return;
 549   }
 550
 551   for(Int_t i=0; i<nDigits; i++){
 552     fDigitsUsed[i]=0 ;
 553   }
 554   Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
 555                      //e.g. first digit in this module, first CPV digit etc.
 556   AliPHOSDigit * digit ;
 557   TArrayI clusterdigitslist(maxNDigits) ;
 558   AliPHOSRecPoint * clu = 0 ;
 559   for(Int_t i=0; i<nDigits; i++){
 560     if(fDigitsUsed[i])
 561       continue ;
 562
 563     digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
 564
 565     clu=0 ;
 566
 567     Int_t index ;
 568
 569     //is this digit so energetic that start cluster?
 570     if (( IsInEmc(digit) &&  digit->GetEnergy() > fEmcClusteringThreshold ) ||
 571         ( IsInCpv(digit) &&  digit->GetEnergy() > fCpvClusteringThreshold ) ) {
 572       Int_t iDigitInCluster = 0 ;
 573
 574       if  ( IsInEmc(digit) ) {
 575         // start a new EMC RecPoint
 576         if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 577           fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
 578
 579         fEMCRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
 580         clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 581         fNumberOfEmcClusters++ ;
 582         clu->AddDigit(*digit, digit->GetEnergy()) ;
 583         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
 584         iDigitInCluster++ ;
 585         fDigitsUsed[i]=kTRUE ;
 586       } else {
 587         // start a new CPV cluster
 588         if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 589           fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
 590
 591         fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
 592         clu =  static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 593         fNumberOfCpvClusters++ ;
 594         clu->AddDigit(*digit, digit->GetEnergy())  ;
 595         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;
 596         iDigitInCluster++ ;
 597         fDigitsUsed[i]=kTRUE ;
 598       } // else
 599
 600       //Now scan remaining digits in list to find neigbours of our seed
 601
 602       AliPHOSDigit * digitN ;
 603       index = 0 ;
 604       while (index < iDigitInCluster){ // scan over digits already in cluster
 605         digit =  static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) )  ;
 606         index++ ;
 607         for(Int_t j=iFirst; j<nDigits; j++){
 608           if(fDigitsUsed[j])
 609             continue ;        //look through remaining digits
 610           digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
 611           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 612           switch (ineb ) {
 613           case -1:   //too early (e.g. previous module), do not look before j at subsequent passes
 614             iFirst=j ;
 615             break ;
 616           case 0 :   // not a neighbour
 617             break ;
 618           case 1 :   // are neighbours
 619             clu->AddDigit(*digitN, digitN->GetEnergy());
 620             clusterdigitslist[iDigitInCluster] = j ;
 621             iDigitInCluster++ ;
 622             fDigitsUsed[j]=kTRUE ;
 623             break ;
 624           case 2 :   // too far from each other
 625             goto endofloop;
 626           } // switch
 627
 628         }
 629
 630         endofloop: ; //scanned all possible neighbours for this digit
 631
 632       } // loop over cluster
 633     } // energy theshold
 634   }
 635
 636 }
 637
 638 //____________________________________________________________________________
 639 void AliPHOSClusterizerv1::MakeUnfolding()
 640 {
 641   // Unfolds clusters using the shape of an ElectroMagnetic shower
 642   // Performs unfolding of all EMC/CPV clusters
 643
 644   // Unfold first EMC clusters
 645   if(fNumberOfEmcClusters > 0){
 646
 647     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 648
 649     Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
 650     Int_t index ;
 651     for(index = 0 ; index < numberofNotUnfolded ; index++){
 652
 653       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
 654       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
 655         break ;
 656
 657       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 658       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 659       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 660       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
 661
 662       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 663         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 664
 665         fEMCRecPoints->Remove(emcRecPoint);
 666         fEMCRecPoints->Compress() ;
 667         index-- ;
 668         fNumberOfEmcClusters -- ;
 669         numberofNotUnfolded-- ;
 670       }
 671       else{
 672         emcRecPoint->SetNExMax(1) ; //Only one local maximum
 673       }
 674
 675       delete[] maxAt ;
 676       delete[] maxAtEnergy ;
 677     }
 678   }
 679   // Unfolding of EMC clusters finished
 680
 681
 682   // Unfold now CPV clusters
 683   if(fNumberOfCpvClusters > 0){
 684
 685     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 686
 687     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
 688     Int_t index ;
 689     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
 690
 691       AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
 692
 693       if(recPoint->GetPHOSMod()> nModulesToUnfold)
 694         break ;
 695
 696       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
 697
 698       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 699       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 700       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 701       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
 702
 703       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 704         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 705         fCPVRecPoints->Remove(emcRecPoint);
 706         fCPVRecPoints->Compress() ;
 707         index-- ;
 708         numberofCpvNotUnfolded-- ;
 709         fNumberOfCpvClusters-- ;
 710       }
 711
 712       delete[] maxAt ;
 713       delete[] maxAtEnergy ;
 714     }
 715   }
 716   //Unfolding of Cpv clusters finished
 717
 718 }
 719
 720 //____________________________________________________________________________
 721 Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
 722 {
 723   // Shape of the shower (see PHOS TDR)
 724   // If you change this function, change also the gradient evaluation in ChiSquare()
 725
 726   //for the moment we neglect dependence on the incident angle.
 727
 728   Double_t r2    = x*x + z*z ;
 729   Double_t r4    = r2*r2 ;
 730   Double_t r295  = TMath::Power(r2, 2.95/2.) ;
 731   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 732   return shape ;
 733 }
 734
 735 //____________________________________________________________________________
 736 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 737                                           Int_t nMax,
 738                                           AliPHOSDigit ** maxAt,
 739                                           Float_t * maxAtEnergy)
 740 {
 741   // Performs the unfolding of a cluster with nMax overlapping showers
 742
 743   Int_t nPar = 3 * nMax ;
 744   Float_t * fitparameters = new Float_t[nPar] ;
 745
 746   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
 747
 748   if( !rv ) {
 749     // Fit failed, return and remove cluster
 750     iniEmc->SetNExMax(-1) ;
 751     delete[] fitparameters ;
 752     return ;
 753   }
 754
 755   // create ufolded rec points and fill them with new energy lists
 756   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
 757   // and later correct this number in acordance with actual energy deposition
 758
 759   Int_t nDigits = iniEmc->GetMultiplicity() ;
 760   Float_t * efit = new Float_t[nDigits] ;
 761   Float_t xDigit=0.,zDigit=0. ;
 762   Float_t xpar=0.,zpar=0.,epar=0.  ;
 763   Int_t relid[4] ;
 764   AliPHOSDigit * digit = 0 ;
 765   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 766
 767   TVector3 vIncid ;
 768
 769   Int_t iparam ;
 770   Int_t iDigit ;
 771   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 772     digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
 773     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 774     fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 775     efit[iDigit] = 0;
 776
 777     iparam = 0 ;
 778     while(iparam < nPar ){
 779       xpar = fitparameters[iparam] ;
 780       zpar = fitparameters[iparam+1] ;
 781       epar = fitparameters[iparam+2] ;
 782       iparam += 3 ;
 783 //      fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
 784 //      efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
 785       efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 786     }
 787   }
 788
 789   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
 790   // so that energy deposited in each cell is distributed betwin new clusters proportionally
 791   // to its contribution to efit
 792
 793   Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
 794   Float_t ratio ;
 795
 796   iparam = 0 ;
 797   while(iparam < nPar ){
 798     xpar = fitparameters[iparam] ;
 799     zpar = fitparameters[iparam+1] ;
 800     epar = fitparameters[iparam+2] ;
 801     iparam += 3 ;
 802 //    fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
 803
 804     AliPHOSEmcRecPoint * emcRP = 0 ;
 805
 806     if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
 807
 808       if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 809         fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
 810
 811       (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
 812       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 813       fNumberOfEmcClusters++ ;
 814       emcRP->SetNExMax((Int_t)nPar/3) ;
 815     }
 816     else{//create new entries in fCPVRecPoints
 817       if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 818         fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
 819
 820       (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
 821       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 822       fNumberOfCpvClusters++ ;
 823     }
 824
 825     Float_t eDigit ;
 826     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 827       digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
 828       fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 829       fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 830 //      ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
 831       ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
 832       eDigit = emcEnergies[iDigit] * ratio ;
 833       emcRP->AddDigit( *digit, eDigit ) ;
 834     }
 835   }
 836
 837   delete[] fitparameters ;
 838   delete[] efit ;
 839
 840 }
 841
 842 //_____________________________________________________________________________
 843 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 844 {
 845   // Calculates the Chi square for the cluster unfolding minimization
 846   // Number of parameters, Gradient, Chi squared, parameters, what to do
 847
 848   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
 849
 850   AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
 851   TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
 852   // A bit buggy way to get an access to the geometry
 853   // To be revised!
 854   AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
 855
 856 //  TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ;  //Vertex position
 857
 858   //  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
 859
 860   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 861
 862   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
 863
 864   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
 865
 866 //  TVector3 vInc ;
 867   fret = 0. ;
 868   Int_t iparam ;
 869
 870   if(iflag == 2)
 871     for(iparam = 0 ; iparam < nPar ; iparam++)
 872       Grad[iparam] = 0 ; // Will evaluate gradient
 873
 874   Double_t efit ;
 875
 876   AliPHOSDigit * digit ;
 877   Int_t iDigit ;
 878
 879   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 880
 881     digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
 882
 883     Int_t relid[4] ;
 884     Float_t xDigit ;
 885     Float_t zDigit ;
 886
 887     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 888
 889     geom->RelPosInModule(relid, xDigit, zDigit) ;
 890
 891      if(iflag == 2){  // calculate gradient
 892        Int_t iParam = 0 ;
 893        efit = 0 ;
 894        while(iParam < nPar ){
 895          Double_t dx = (xDigit - x[iParam]) ;
 896          iParam++ ;
 897          Double_t dz = (zDigit - x[iParam]) ;
 898          iParam++ ;
 899 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
 900 //         efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
 901          efit += x[iParam] * ShowerShape(dx,dz) ;
 902          iParam++ ;
 903        }
 904        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
 905        iParam = 0 ;
 906        while(iParam < nPar ){
 907          Double_t xpar = x[iParam] ;
 908          Double_t zpar = x[iParam+1] ;
 909          Double_t epar = x[iParam+2] ;
 910 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
 911          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
 912 //         Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
 913          Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 914 //DP: No incident angle dependence in gradient yet!!!!!!
 915          Double_t r4 = dr*dr*dr*dr ;
 916          Double_t r295 = TMath::Power(dr,2.95) ;
 917          Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
 918                                          0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
 919
 920          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
 921          iParam++ ;
 922          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
 923          iParam++ ;
 924          Grad[iParam] += shape ;                                  // Derivative over energy
 925          iParam++ ;
 926        }
 927      }
 928      efit = 0;
 929      iparam = 0 ;
 930
 931      while(iparam < nPar ){
 932        Double_t xpar = x[iparam] ;
 933        Double_t zpar = x[iparam+1] ;
 934        Double_t epar = x[iparam+2] ;
 935        iparam += 3 ;
 936 //       fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
 937 //       efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
 938        efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 939      }
 940
 941      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
 942      // Here we assume, that sigma = sqrt(E)
 943   }
 944
 945 }
 946
 947 //____________________________________________________________________________
 948 void AliPHOSClusterizerv1::Print(const Option_t *)const
 949 {
 950   // Print clusterizer parameters
 951
 952   TString message ;
 953   TString taskName(GetName()) ;
 954   taskName.ReplaceAll(Version(), "") ;
 955
 956   if( strcmp(GetName(), "") !=0 ) {
 957     // Print parameters
 958     message  = "\n--------------- %s %s -----------\n" ;
 959     message += "Clusterizing digits from the file: %s\n" ;
 960     message += "                           Branch: %s\n" ;
 961     message += "                       EMC Clustering threshold = %f\n" ;
 962     message += "                       EMC Local Maximum cut    = %f\n" ;
 963     message += "                       EMC Logarothmic weight   = %f\n" ;
 964     message += "                       CPV Clustering threshold = %f\n" ;
 965     message += "                       CPV Local Maximum cut    = %f\n" ;
 966     message += "                       CPV Logarothmic weight   = %f\n" ;
 967     if(fToUnfold)
 968       message += " Unfolding on\n" ;
 969     else
 970       message += " Unfolding off\n" ;
 971
 972     message += "------------------------------------------------------------------" ;
 973   }
 974   else
 975     message = " AliPHOSClusterizerv1 not initialized " ;
 976
 977   AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
 978        taskName.Data(),
 979        GetTitle(),
 980        taskName.Data(),
 981        GetName(),
 982        fEmcClusteringThreshold,
 983        fEmcLocMaxCut,
 984        fW0,
 985        fCpvClusteringThreshold,
 986        fCpvLocMaxCut,
 987        fW0CPV )) ;
 988 }
 989 //____________________________________________________________________________
 990 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
 991 {
 992   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
 993
 994   AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
 995                fEMCRecPoints->GetEntriesFast(),
 996                fCPVRecPoints->GetEntriesFast() ))  ;
 997
 998   if(strstr(option,"all")) {
 999     printf("\n EMC clusters \n") ;
1000     printf("Index    Ene(MeV) Multi Module     X    Y   Z    Lambdas_1  Lambda_2  # of prim  Primaries list\n") ;
1001     Int_t index ;
1002     for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1003       AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1004       TVector3  locpos;
1005       rp->GetLocalPosition(locpos);
1006       Float_t lambda[2];
1007       rp->GetElipsAxis(lambda);
1008       Int_t * primaries;
1009       Int_t nprimaries;
1010       primaries = rp->GetPrimaries(nprimaries);
1011       printf("\n%6d  %8.2f  %3d     %2d     %4.1f   %4.1f   %4.1f   %4f  %4f    %2d     : ",
1012               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1013               locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1014
1015       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1016         printf("%d ", primaries[iprimary] ) ;
1017       }
1018       printf("\n") ;
1019     }
1020
1021     //Now plot CPV recPoints
1022     printf("\n CPV clusters \n") ;
1023     printf("Index    Ene(MeV) Module     X     Y    Z  \n") ;
1024     for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1025       AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1026
1027       TVector3  locpos;
1028       rp->GetLocalPosition(locpos);
1029
1030       printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f \n",
1031              rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1032              locpos.X(), locpos.Y(), locpos.Z()) ;
1033     }
1034   }
1035 }
1036
1037
1038 //____________________________________________________________________________
1039 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1040 {
1041   //For each EMC rec. point set the distance to the nearest bad crystal.
1042   //Author: Boris Polichtchouk
1043
1044   if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1045   AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
1046
1047   Int_t badIds[8000];
1048   fgCalibData->EmcBadChannelIds(badIds);
1049
1050   AliPHOSEmcRecPoint* rp;
1051
1052   TMatrixF gmat;
1053   TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1054   TVector3 gposBadChannel; // global position of bad crystal
1055   TVector3 dR;
1056
1057   Float_t dist,minDist;
1058   Int_t relid[4]={0,0,0,0} ;
1059   TVector3 lpos ;
1060   for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1061     rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1062     //evaluate distance to border
1063     relid[0]=rp->GetPHOSMod() ;
1064     relid[2]=1 ;
1065     relid[3]=1 ;
1066     Float_t xcorner,zcorner;
1067     fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1068     rp->GetLocalPosition(lpos) ;
1069     minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
1070     for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1071       fGeom->AbsToRelNumbering(badIds[iBad],relid)  ;
1072       if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1073         continue ;                   //bad channels can be in the module which does not exist in simulations.
1074       rp->GetGlobalPosition(gposRecPoint,gmat);
1075       fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1076       AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1077                       gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1078                       gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1079       dR = gposBadChannel-gposRecPoint;
1080       dist = dR.Mag();
1081       if(dist<minDist) minDist = dist;
1082     }
1083
1084     rp->SetDistanceToBadCrystal(minDist);
1085   }
1086
1087 }