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