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   InitParameters() ;
 206   fDefaultInit = kTRUE ;
 207 }
 208
 209 //____________________________________________________________________________
 210 AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
 211   AliPHOSClusterizer(geom),
 212   fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
 213   fWrite(0),
 214   fNumberOfEmcClusters(0),    fNumberOfCpvClusters(0),
 215   fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
 216   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
 217   fW0CPV(0),                  fEmcTimeGate(0),          fEcoreRadius(0)
 218 {
 219   // ctor with the indication of the file where header Tree and digits Tree are stored
 220
 221   InitParameters() ;
 222   Init() ;
 223   fDefaultInit = kFALSE ;
 224 }
 225
 226 //____________________________________________________________________________
 227   AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
 228 {
 229   // dtor
 230
 231 }
 232 //____________________________________________________________________________
 233 void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
 234 {
 235   // Steering method to perform clusterization for one event
 236   // The input is the tree with digits.
 237   // The output is the tree with clusters.
 238
 239   if(strstr(option,"tim"))
 240     gBenchmark->Start("PHOSClusterizer");
 241
 242   if(strstr(option,"print")) {
 243     Print() ;
 244     return ;
 245   }
 246
 247   MakeClusters() ;
 248
 249   AliDebug(2,Form(" ---- Printing clusters (%d)\n",
 250                   fEMCRecPoints->GetEntries()));
 251   if(AliLog::GetGlobalDebugLevel()>1)
 252     fEMCRecPoints->Print();
 253
 254   if(fToUnfold)
 255     MakeUnfolding();
 256
 257   WriteRecPoints();
 258
 259   if(strstr(option,"deb"))
 260     PrintRecPoints(option) ;
 261
 262   if(strstr(option,"tim")){
 263     gBenchmark->Stop("PHOSClusterizer");
 264     AliInfo(Form("took %f seconds for Clusterizing\n",
 265                  gBenchmark->GetCpuTime("PHOSClusterizer")));
 266   }
 267   fEMCRecPoints->Delete();
 268   fCPVRecPoints->Delete();
 269 }
 270
 271 //____________________________________________________________________________
 272 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
 273                                     Int_t nPar, Float_t * fitparameters) const
 274 {
 275   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 276   // The initial values for fitting procedure are set equal to the positions of local maxima.
 277   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 278
 279
 280   gMinuit->mncler();                     // Reset Minuit's list of paramters
 281   gMinuit->SetPrintLevel(-1) ;           // No Printout
 282   gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
 283                                          // To set the address of the minimization function
 284
 285   TList * toMinuit = new TList();
 286   toMinuit->AddAt(emcRP,0) ;
 287   toMinuit->AddAt(fDigitsArr,1) ;
 288   toMinuit->AddAt(fGeom,2) ;
 289
 290   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 291
 292   // filling initial values for fit parameters
 293   AliPHOSDigit * digit ;
 294
 295   Int_t ierflg  = 0;
 296   Int_t index   = 0 ;
 297   Int_t nDigits = (Int_t) nPar / 3 ;
 298
 299   Int_t iDigit ;
 300
 301   for(iDigit = 0; iDigit < nDigits; iDigit++){
 302     digit = maxAt[iDigit];
 303
 304     Int_t relid[4] ;
 305     Float_t x = 0.;
 306     Float_t z = 0.;
 307     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 308     fGeom->RelPosInModule(relid, x, z) ;
 309
 310     Float_t energy = maxAtEnergy[iDigit] ;
 311
 312     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 313     index++ ;
 314     if(ierflg != 0){
 315       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
 316       return kFALSE;
 317     }
 318     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 319     index++ ;
 320     if(ierflg != 0){
 321        Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
 322       return kFALSE;
 323     }
 324     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 325     index++ ;
 326     if(ierflg != 0){
 327       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
 328       return kFALSE;
 329     }
 330   }
 331
 332   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 333                       //  depends on it.
 334   Double_t p1 = 1.0 ;
 335   Double_t p2 = 0.0 ;
 336
 337   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 338   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 339   gMinuit->SetMaxIterations(5);
 340   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 341
 342   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 343
 344   if(ierflg == 4){  // Minimum not found
 345     Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
 346     return kFALSE ;
 347   }
 348   for(index = 0; index < nPar; index++){
 349     Double_t err ;
 350     Double_t val ;
 351     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 352     fitparameters[index] = val ;
 353    }
 354
 355   delete toMinuit ;
 356   return kTRUE;
 357
 358 }
 359
 360
 361 //____________________________________________________________________________
 362 void AliPHOSClusterizerv1::Init()
 363 {
 364   // Make all memory allocations which can not be done in default constructor.
 365   // Attach the Clusterizer task to the list of PHOS tasks
 366
 367   fEmcCrystals = fGeom->GetNModules() *  fGeom->GetNCristalsInModule() ;
 368
 369   if(!gMinuit)
 370     gMinuit = new TMinuit(100);
 371
 372   if (!fgCalibData)
 373     fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
 374   if (fgCalibData->GetCalibDataEmc() == 0)
 375     AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
 376
 377 }
 378
 379 //____________________________________________________________________________
 380 void AliPHOSClusterizerv1::InitParameters()
 381 {
 382
 383   fNumberOfCpvClusters     = 0 ;
 384   fNumberOfEmcClusters     = 0 ;
 385
 386   const AliPHOSRecoParam* parEmc = AliPHOSReconstructor::GetRecoParamEmc();
 387   if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!");
 388
 389   const AliPHOSRecoParam* parCpv = AliPHOSReconstructor::GetRecoParamCpv();
 390   if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!");
 391
 392   fCpvClusteringThreshold  = parCpv->GetClusteringThreshold();
 393   fEmcClusteringThreshold  = parEmc->GetClusteringThreshold();
 394
 395   fEmcLocMaxCut            = parEmc->GetLocalMaxCut();
 396   fCpvLocMaxCut            = parCpv->GetLocalMaxCut();
 397
 398   fW0                      = parEmc->GetLogWeight();
 399   fW0CPV                   = parCpv->GetLogWeight();
 400
 401   fEmcTimeGate             = 1.e-6 ;
 402   fEcoreRadius             = parEmc->GetEcoreRadius();
 403
 404   fToUnfold                = parEmc->ToUnfold() ;
 405
 406   fWrite                   = kTRUE ;
 407 }
 408
 409 //____________________________________________________________________________
 410 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
 411 {
 412   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 413   //                                       = 1 are neighbour
 414   //                                       = 2 are not neighbour but do not continue searching
 415   //                                       =-1 are not neighbour, continue searching, but do not look before d2 next time
 416   // neighbours are defined as digits having at least a common vertex
 417   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 418   //                                      which is compared to a digit (d2)  not yet in a cluster
 419
 420   Int_t relid1[4] ;
 421   fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
 422
 423   Int_t relid2[4] ;
 424   fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
 425
 426   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
 427     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 428     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 429
 430     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){   //At least common vertex
 431       //    if (( relid1[2]==relid2[2] && coldiff <= 1 )  || ( relid1[3]==relid2[3] &&  rowdiff <= 1 )){ //common side
 432       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
 433       return 1 ;
 434     }
 435     else {
 436       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
 437         return 2; //  Difference in row numbers is too large to look further
 438     }
 439     return 0 ;
 440
 441   }
 442   else {
 443     if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
 444       return -1 ;
 445     if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
 446       return -1 ;
 447
 448     return 2 ;
 449
 450   }
 451
 452   return 0 ;
 453 }
 454 //____________________________________________________________________________
 455 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
 456 {
 457   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
 458
 459   Bool_t rv = kFALSE ;
 460
 461   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 462
 463   if(digit->GetId() <= nEMC )   rv = kTRUE;
 464
 465   return rv ;
 466 }
 467
 468 //____________________________________________________________________________
 469 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 470 {
 471   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
 472
 473   Bool_t rv = kFALSE ;
 474
 475   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 476
 477   if(digit->GetId() > nEMC )   rv = kTRUE;
 478
 479   return rv ;
 480 }
 481
 482 //____________________________________________________________________________
 483 void AliPHOSClusterizerv1::WriteRecPoints()
 484 {
 485
 486   // Creates new branches with given title
 487   // fills and writes into TreeR.
 488
 489   Int_t index ;
 490   //Evaluate position, dispersion and other RecPoint properties..
 491   Int_t nEmc = fEMCRecPoints->GetEntriesFast();
 492   Float_t emcMinE= AliPHOSReconstructor::GetRecoParamEmc()->GetMinE(); //Minimal digit energy
 493   TVector3 fakeVtx(0.,0.,0.) ;
 494   for(index = 0; index < nEmc; index++){
 495     AliPHOSEmcRecPoint * rp =
 496       dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
 497     rp->Purify(emcMinE) ;
 498     if(rp->GetMultiplicity()==0){
 499       fEMCRecPoints->RemoveAt(index) ;
 500       delete rp ;
 501       continue;
 502     }
 503
 504     // No vertex is available now, calculate corrections in PID
 505     rp->EvalAll(fDigitsArr) ;
 506     rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
 507     rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
 508     rp->EvalLocal2TrackingCSTransform();
 509   }
 510   fEMCRecPoints->Compress() ;
 511   fEMCRecPoints->Sort() ;
 512   //  fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
 513   for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
 514     dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
 515   }
 516
 517   //For each rec.point set the distance to the nearest bad crystal (BVP)
 518   SetDistancesToBadChannels();
 519
 520   //Now the same for CPV
 521   for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
 522     AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
 523     rp->EvalAll(fDigitsArr) ;
 524     rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
 525     rp->EvalLocal2TrackingCSTransform();
 526   }
 527   fCPVRecPoints->Sort() ;
 528
 529   for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
 530     dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
 531
 532   fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
 533
 534   if(fWrite){ //We write TreeR
 535     fTreeR->Fill();
 536   }
 537 }
 538
 539 //____________________________________________________________________________
 540 void AliPHOSClusterizerv1::MakeClusters()
 541 {
 542   // Steering method to construct the clusters stored in a list of Reconstructed Points
 543   // A cluster is defined as a list of neighbour digits
 544
 545   fNumberOfCpvClusters     = 0 ;
 546   fNumberOfEmcClusters     = 0 ;
 547
 548   //Mark all digits as unused yet
 549   Int_t nDigits=fDigitsArr->GetEntriesFast() ;
 550   for(Int_t i=0; i<nDigits; i++){
 551     fDigitsUsed[i]=0 ;
 552   }
 553   Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
 554                      //e.g. first digit in this module, first CPV digit etc.
 555   AliPHOSDigit * digit ;
 556   TArrayI clusterdigitslist(1500) ;
 557   AliPHOSRecPoint * clu = 0 ;
 558   for(Int_t i=0; i<nDigits; i++){
 559     if(fDigitsUsed[i])
 560       continue ;
 561
 562     digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
 563
 564     clu=0 ;
 565
 566     Int_t index ;
 567
 568     //is this digit so energetic that start cluster?
 569     if (( IsInEmc(digit) &&  digit->GetEnergy() > fEmcClusteringThreshold ) ||
 570         ( IsInCpv(digit) &&  digit->GetEnergy() > fCpvClusteringThreshold ) ) {
 571       Int_t iDigitInCluster = 0 ;
 572
 573       if  ( IsInEmc(digit) ) {
 574         // start a new EMC RecPoint
 575         if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 576           fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
 577
 578         fEMCRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
 579         clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 580         fNumberOfEmcClusters++ ;
 581         clu->AddDigit(*digit, digit->GetEnergy()) ;
 582         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
 583         iDigitInCluster++ ;
 584         fDigitsUsed[i]=kTRUE ;
 585       } else {
 586         // start a new CPV cluster
 587         if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 588           fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
 589
 590         fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
 591         clu =  static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 592         fNumberOfCpvClusters++ ;
 593         clu->AddDigit(*digit, digit->GetEnergy())  ;
 594         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;
 595         iDigitInCluster++ ;
 596         fDigitsUsed[i]=kTRUE ;
 597       } // else
 598
 599       //Now scan remaining digits in list to find neigbours of our seed
 600
 601       AliPHOSDigit * digitN ;
 602       index = 0 ;
 603       while (index < iDigitInCluster){ // scan over digits already in cluster
 604         digit =  static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) )  ;
 605         index++ ;
 606         for(Int_t j=iFirst; j<nDigits; j++){
 607           if(fDigitsUsed[j])
 608             continue ;        //look through remaining digits
 609           digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
 610           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 611           switch (ineb ) {
 612           case -1:   //too early (e.g. previous module), do not look before j at subsequent passes
 613             iFirst=j ;
 614             break ;
 615           case 0 :   // not a neighbour
 616             break ;
 617           case 1 :   // are neighbours
 618             clu->AddDigit(*digitN, digitN->GetEnergy());
 619             clusterdigitslist[iDigitInCluster] = j ;
 620             iDigitInCluster++ ;
 621             fDigitsUsed[j]=kTRUE ;
 622             break ;
 623           case 2 :   // too far from each other
 624             goto endofloop;
 625           } // switch
 626
 627         }
 628
 629         endofloop: ; //scanned all possible neighbours for this digit
 630
 631       } // loop over cluster
 632     } // energy theshold
 633   }
 634
 635 }
 636
 637 //____________________________________________________________________________
 638 void AliPHOSClusterizerv1::MakeUnfolding()
 639 {
 640   // Unfolds clusters using the shape of an ElectroMagnetic shower
 641   // Performs unfolding of all EMC/CPV clusters
 642
 643   // Unfold first EMC clusters
 644   if(fNumberOfEmcClusters > 0){
 645
 646     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 647
 648     Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
 649     Int_t index ;
 650     for(index = 0 ; index < numberofNotUnfolded ; index++){
 651
 652       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
 653       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
 654         break ;
 655
 656       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 657       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 658       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 659       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
 660
 661       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 662         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 663
 664         fEMCRecPoints->Remove(emcRecPoint);
 665         fEMCRecPoints->Compress() ;
 666         index-- ;
 667         fNumberOfEmcClusters -- ;
 668         numberofNotUnfolded-- ;
 669       }
 670       else{
 671         emcRecPoint->SetNExMax(1) ; //Only one local maximum
 672       }
 673
 674       delete[] maxAt ;
 675       delete[] maxAtEnergy ;
 676     }
 677   }
 678   // Unfolding of EMC clusters finished
 679
 680
 681   // Unfold now CPV clusters
 682   if(fNumberOfCpvClusters > 0){
 683
 684     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 685
 686     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
 687     Int_t index ;
 688     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
 689
 690       AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
 691
 692       if(recPoint->GetPHOSMod()> nModulesToUnfold)
 693         break ;
 694
 695       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
 696
 697       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 698       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 699       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 700       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
 701
 702       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 703         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 704         fCPVRecPoints->Remove(emcRecPoint);
 705         fCPVRecPoints->Compress() ;
 706         index-- ;
 707         numberofCpvNotUnfolded-- ;
 708         fNumberOfCpvClusters-- ;
 709       }
 710
 711       delete[] maxAt ;
 712       delete[] maxAtEnergy ;
 713     }
 714   }
 715   //Unfolding of Cpv clusters finished
 716
 717 }
 718
 719 //____________________________________________________________________________
 720 Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
 721 {
 722   // Shape of the shower (see PHOS TDR)
 723   // If you change this function, change also the gradient evaluation in ChiSquare()
 724
 725   //for the moment we neglect dependence on the incident angle.
 726
 727   Double_t r2    = x*x + z*z ;
 728   Double_t r4    = r2*r2 ;
 729   Double_t r295  = TMath::Power(r2, 2.95/2.) ;
 730   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 731   return shape ;
 732 }
 733
 734 //____________________________________________________________________________
 735 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 736                                           Int_t nMax,
 737                                           AliPHOSDigit ** maxAt,
 738                                           Float_t * maxAtEnergy)
 739 {
 740   // Performs the unfolding of a cluster with nMax overlapping showers
 741
 742   Int_t nPar = 3 * nMax ;
 743   Float_t * fitparameters = new Float_t[nPar] ;
 744
 745   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
 746
 747   if( !rv ) {
 748     // Fit failed, return and remove cluster
 749     iniEmc->SetNExMax(-1) ;
 750     delete[] fitparameters ;
 751     return ;
 752   }
 753
 754   // create ufolded rec points and fill them with new energy lists
 755   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
 756   // and later correct this number in acordance with actual energy deposition
 757
 758   Int_t nDigits = iniEmc->GetMultiplicity() ;
 759   Float_t * efit = new Float_t[nDigits] ;
 760   Float_t xDigit=0.,zDigit=0. ;
 761   Float_t xpar=0.,zpar=0.,epar=0.  ;
 762   Int_t relid[4] ;
 763   AliPHOSDigit * digit = 0 ;
 764   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 765
 766   TVector3 vIncid ;
 767
 768   Int_t iparam ;
 769   Int_t iDigit ;
 770   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 771     digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
 772     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 773     fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 774     efit[iDigit] = 0;
 775
 776     iparam = 0 ;
 777     while(iparam < nPar ){
 778       xpar = fitparameters[iparam] ;
 779       zpar = fitparameters[iparam+1] ;
 780       epar = fitparameters[iparam+2] ;
 781       iparam += 3 ;
 782 //      fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
 783 //      efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
 784       efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 785     }
 786   }
 787
 788   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
 789   // so that energy deposited in each cell is distributed betwin new clusters proportionally
 790   // to its contribution to efit
 791
 792   Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
 793   Float_t ratio ;
 794
 795   iparam = 0 ;
 796   while(iparam < nPar ){
 797     xpar = fitparameters[iparam] ;
 798     zpar = fitparameters[iparam+1] ;
 799     epar = fitparameters[iparam+2] ;
 800     iparam += 3 ;
 801 //    fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
 802
 803     AliPHOSEmcRecPoint * emcRP = 0 ;
 804
 805     if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
 806
 807       if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 808         fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
 809
 810       (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
 811       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 812       fNumberOfEmcClusters++ ;
 813       emcRP->SetNExMax((Int_t)nPar/3) ;
 814     }
 815     else{//create new entries in fCPVRecPoints
 816       if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 817         fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
 818
 819       (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
 820       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 821       fNumberOfCpvClusters++ ;
 822     }
 823
 824     Float_t eDigit ;
 825     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 826       digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
 827       fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 828       fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 829 //      ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
 830       ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
 831       eDigit = emcEnergies[iDigit] * ratio ;
 832       emcRP->AddDigit( *digit, eDigit ) ;
 833     }
 834   }
 835
 836   delete[] fitparameters ;
 837   delete[] efit ;
 838
 839 }
 840
 841 //_____________________________________________________________________________
 842 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 843 {
 844   // Calculates the Chi square for the cluster unfolding minimization
 845   // Number of parameters, Gradient, Chi squared, parameters, what to do
 846
 847   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
 848
 849   AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
 850   TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
 851   // A bit buggy way to get an access to the geometry
 852   // To be revised!
 853   AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
 854
 855 //  TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ;  //Vertex position
 856
 857   //  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
 858
 859   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 860
 861   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
 862
 863   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
 864
 865 //  TVector3 vInc ;
 866   fret = 0. ;
 867   Int_t iparam ;
 868
 869   if(iflag == 2)
 870     for(iparam = 0 ; iparam < nPar ; iparam++)
 871       Grad[iparam] = 0 ; // Will evaluate gradient
 872
 873   Double_t efit ;
 874
 875   AliPHOSDigit * digit ;
 876   Int_t iDigit ;
 877
 878   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 879
 880     digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
 881
 882     Int_t relid[4] ;
 883     Float_t xDigit ;
 884     Float_t zDigit ;
 885
 886     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 887
 888     geom->RelPosInModule(relid, xDigit, zDigit) ;
 889
 890      if(iflag == 2){  // calculate gradient
 891        Int_t iParam = 0 ;
 892        efit = 0 ;
 893        while(iParam < nPar ){
 894          Double_t dx = (xDigit - x[iParam]) ;
 895          iParam++ ;
 896          Double_t dz = (zDigit - x[iParam]) ;
 897          iParam++ ;
 898 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
 899 //         efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
 900          efit += x[iParam] * ShowerShape(dx,dz) ;
 901          iParam++ ;
 902        }
 903        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
 904        iParam = 0 ;
 905        while(iParam < nPar ){
 906          Double_t xpar = x[iParam] ;
 907          Double_t zpar = x[iParam+1] ;
 908          Double_t epar = x[iParam+2] ;
 909 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
 910          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
 911 //         Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
 912          Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 913 //DP: No incident angle dependence in gradient yet!!!!!!
 914          Double_t r4 = dr*dr*dr*dr ;
 915          Double_t r295 = TMath::Power(dr,2.95) ;
 916          Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
 917                                          0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
 918
 919          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
 920          iParam++ ;
 921          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
 922          iParam++ ;
 923          Grad[iParam] += shape ;                                  // Derivative over energy
 924          iParam++ ;
 925        }
 926      }
 927      efit = 0;
 928      iparam = 0 ;
 929
 930      while(iparam < nPar ){
 931        Double_t xpar = x[iparam] ;
 932        Double_t zpar = x[iparam+1] ;
 933        Double_t epar = x[iparam+2] ;
 934        iparam += 3 ;
 935 //       fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
 936 //       efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
 937        efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 938      }
 939
 940      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
 941      // Here we assume, that sigma = sqrt(E)
 942   }
 943
 944 }
 945
 946 //____________________________________________________________________________
 947 void AliPHOSClusterizerv1::Print(const Option_t *)const
 948 {
 949   // Print clusterizer parameters
 950
 951   TString message ;
 952   TString taskName(GetName()) ;
 953   taskName.ReplaceAll(Version(), "") ;
 954
 955   if( strcmp(GetName(), "") !=0 ) {
 956     // Print parameters
 957     message  = "\n--------------- %s %s -----------\n" ;
 958     message += "Clusterizing digits from the file: %s\n" ;
 959     message += "                           Branch: %s\n" ;
 960     message += "                       EMC Clustering threshold = %f\n" ;
 961     message += "                       EMC Local Maximum cut    = %f\n" ;
 962     message += "                       EMC Logarothmic weight   = %f\n" ;
 963     message += "                       CPV Clustering threshold = %f\n" ;
 964     message += "                       CPV Local Maximum cut    = %f\n" ;
 965     message += "                       CPV Logarothmic weight   = %f\n" ;
 966     if(fToUnfold)
 967       message += " Unfolding on\n" ;
 968     else
 969       message += " Unfolding off\n" ;
 970
 971     message += "------------------------------------------------------------------" ;
 972   }
 973   else
 974     message = " AliPHOSClusterizerv1 not initialized " ;
 975
 976   AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
 977        taskName.Data(),
 978        GetTitle(),
 979        taskName.Data(),
 980        GetName(),
 981        fEmcClusteringThreshold,
 982        fEmcLocMaxCut,
 983        fW0,
 984        fCpvClusteringThreshold,
 985        fCpvLocMaxCut,
 986        fW0CPV )) ;
 987 }
 988 //____________________________________________________________________________
 989 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
 990 {
 991   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
 992
 993   AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
 994                fEMCRecPoints->GetEntriesFast(),
 995                fCPVRecPoints->GetEntriesFast() ))  ;
 996
 997   if(strstr(option,"all")) {
 998     printf("\n EMC clusters \n") ;
 999     printf("Index    Ene(MeV) Multi Module     X    Y   Z    Lambdas_1  Lambda_2  # of prim  Primaries list\n") ;
1000     Int_t index ;
1001     for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1002       AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1003       TVector3  locpos;
1004       rp->GetLocalPosition(locpos);
1005       Float_t lambda[2];
1006       rp->GetElipsAxis(lambda);
1007       Int_t * primaries;
1008       Int_t nprimaries;
1009       primaries = rp->GetPrimaries(nprimaries);
1010       printf("\n%6d  %8.2f  %3d     %2d     %4.1f   %4.1f   %4.1f   %4f  %4f    %2d     : ",
1011               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1012               locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1013
1014       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1015         printf("%d ", primaries[iprimary] ) ;
1016       }
1017       printf("\n") ;
1018     }
1019
1020     //Now plot CPV recPoints
1021     printf("\n CPV clusters \n") ;
1022     printf("Index    Ene(MeV) Module     X     Y    Z  \n") ;
1023     for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1024       AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1025
1026       TVector3  locpos;
1027       rp->GetLocalPosition(locpos);
1028
1029       printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f \n",
1030              rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1031              locpos.X(), locpos.Y(), locpos.Z()) ;
1032     }
1033   }
1034 }
1035
1036
1037 //____________________________________________________________________________
1038 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1039 {
1040   //For each EMC rec. point set the distance to the nearest bad crystal.
1041   //Author: Boris Polichtchouk
1042
1043   if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1044   AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
1045
1046   Int_t badIds[8000];
1047   fgCalibData->EmcBadChannelIds(badIds);
1048
1049   AliPHOSEmcRecPoint* rp;
1050
1051   TMatrixF gmat;
1052   TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1053   TVector3 gposBadChannel; // global position of bad crystal
1054   TVector3 dR;
1055
1056   Float_t dist,minDist;
1057   Int_t relid[4]={0,0,0,0} ;
1058   TVector3 lpos ;
1059   for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1060     rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1061     //evaluate distance to border
1062     relid[0]=rp->GetPHOSMod() ;
1063     relid[2]=1 ;
1064     relid[3]=1 ;
1065     Float_t xcorner,zcorner;
1066     fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1067     rp->GetLocalPosition(lpos) ;
1068     minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
1069     for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1070       fGeom->AbsToRelNumbering(badIds[iBad],relid)  ;
1071       if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1072         continue ;                   //bad channels can be in the module which does not exist in simulations.
1073       rp->GetGlobalPosition(gposRecPoint,gmat);
1074       fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1075       AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1076                       gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1077                       gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1078       dR = gposBadChannel-gposRecPoint;
1079       dist = dR.Mag();
1080       if(dist<minDist) minDist = dist;
1081     }
1082
1083     rp->SetDistanceToBadCrystal(minDist);
1084   }
1085
1086 }