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