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