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