PHOS/AliPHOSClusterizerv1.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /* History of cvs commits:
  19  *
  20  * $Log: AliPHOSClusterizerv1.cxx,v $
  21  * Revision 1.118  2007/12/11 21:23:26  kharlov
  22  * Added possibility to swith off unfolding
  23  *
  24  * Revision 1.117  2007/10/18 08:42:05  kharlov
  25  * Bad channels cleaned before clusterization
  26  *
  27  * Revision 1.116  2007/10/01 20:24:08  kharlov
  28  * Memory leaks fixed
  29  *
  30  * Revision 1.115  2007/09/26 14:22:17  cvetan
  31  * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
  32  *
  33  * Revision 1.114  2007/09/06 16:06:44  kharlov
  34  * Absence of sorting results in loose of all unfolded clusters
  35  *
  36  * Revision 1.113  2007/08/28 12:55:07  policheh
  37  * Loaders removed from the reconstruction code (C.Cheshkov)
  38  *
  39  * Revision 1.112  2007/08/22 09:20:50  hristov
  40  * Updated QA classes (Yves)
  41  *
  42  * Revision 1.111  2007/08/08 12:11:28  kharlov
  43  * Protection against uninitialized fQADM
  44  *
  45  * Revision 1.110  2007/08/07 14:16:00  kharlov
  46  * Quality assurance added (Yves Schutz)
  47  *
  48  * Revision 1.109  2007/07/24 17:20:35  policheh
  49  * Usage of RecoParam objects instead of hardcoded parameters in reconstruction.
  50  * (See $ALICE_ROOT/PHOS/macros/BeamTest2006/RawReconstruction.C).
  51  *
  52  * Revision 1.108  2007/06/18 07:00:51  kharlov
  53  * Bug fix for attempt to use AliPHOSEmcRecPoint after its deletion
  54  *
  55  * Revision 1.107  2007/05/25 14:12:26  policheh
  56  * Local to tracking CS transformation added for CPV rec. points
  57  *
  58  * Revision 1.106  2007/05/24 13:01:22  policheh
  59  * Local to tracking CS transformation invoked for each EMC rec.point
  60  *
  61  * Revision 1.105  2007/05/02 13:41:22  kharlov
  62  * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
  63  *
  64  * Revision 1.104  2007/04/27 16:55:53  kharlov
  65  * Calibration stops if PHOS CDB objects do not exist
  66  *
  67  * Revision 1.103  2007/04/11 11:55:45  policheh
  68  * SetDistancesToBadChannels() added.
  69  *
  70  * Revision 1.102  2007/03/28 19:18:15  kharlov
  71  * RecPoints recalculation in TSM removed
  72  *
  73  * Revision 1.101  2007/03/06 06:51:27  kharlov
  74  * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
  75  *
  76  * Revision 1.100  2007/01/10 11:07:26  kharlov
  77  * Raw digits writing to file (B.Polichtchouk)
  78  *
  79  * Revision 1.99  2006/11/07 16:49:51  kharlov
  80  * Corrections for next event switch in case of raw data (B.Polichtchouk)
  81  *
  82  * Revision 1.98  2006/10/27 17:14:27  kharlov
  83  * Introduce AliDebug and AliLog (B.Polichtchouk)
  84  *
  85  * Revision 1.97  2006/08/29 11:41:19  kharlov
  86  * Missing implementation of ctors and = operator are added
  87  *
  88  * Revision 1.96  2006/08/25 16:56:30  kharlov
  89  * Compliance with Effective C++
  90  *
  91  * Revision 1.95  2006/08/11 12:36:26  cvetan
  92  * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
  93  *
  94  * Revision 1.94  2006/08/07 12:27:49  hristov
  95  * Removing obsolete code which affected the event numbering scheme
  96  *
  97  * Revision 1.93  2006/08/01 12:20:17  cvetan
  98  * 1. Adding a possibility to read and reconstruct an old rcu formatted raw data. This is controlled by an option of AliReconstruction and AliPHOSReconstructor. 2. In case of raw data processing (without galice.root) create the default AliPHOSGeometry object. Most likely this should be moved to the CDB
  99  *
 100  * Revision 1.92  2006/04/29 20:26:46  hristov
 101  * Separate EMC and CPV calibration (Yu.Kharlov)
 102  *
 103  * Revision 1.91  2006/04/22 10:30:17  hristov
 104  * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
 105  *
 106  * Revision 1.90  2006/04/11 15:22:59  hristov
 107  * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
 108  *
 109  * Revision 1.89  2006/03/13 14:05:42  kharlov
 110  * Calibration objects for EMC and CPV
 111  *
 112  * Revision 1.88  2006/01/11 08:54:52  hristov
 113  * Additional protection in case no calibration entry was found
 114  *
 115  * Revision 1.87  2005/11/22 08:46:43  kharlov
 116  * Updated to new CDB (Boris Polichtchouk)
 117  *
 118  * Revision 1.86  2005/11/14 21:52:43  hristov
 119  * Coding conventions
 120  *
 121  * Revision 1.85  2005/09/27 16:08:08  hristov
 122  * New version of CDB storage framework (A.Colla)
 123  *
 124  * Revision 1.84  2005/09/21 10:02:47  kharlov
 125  * Reading calibration from CDB (Boris Polichtchouk)
 126  *
 127  * Revision 1.82  2005/09/02 15:43:13  kharlov
 128  * Add comments in GetCalibrationParameters and Calibrate
 129  *
 130  * Revision 1.81  2005/09/02 14:32:07  kharlov
 131  * Calibration of raw data
 132  *
 133  * Revision 1.80  2005/08/24 15:31:36  kharlov
 134  * Setting raw digits flag
 135  *
 136  * Revision 1.79  2005/07/25 15:53:53  kharlov
 137  * Read raw data
 138  *
 139  * Revision 1.78  2005/05/28 14:19:04  schutz
 140  * Compilation warnings fixed by T.P.
 141  *
 142  */
 143
 144 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
 145 //////////////////////////////////////////////////////////////////////////////
 146 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and
 147 //  unfolds the clusters having several local maxima.
 148 //  Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
 149 //  PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
 150 //  parameters including input digits branch title, thresholds etc.)
 151 //  This TTask is normally called from Reconstructioner, but can as well be used in
 152 //  standalone mode.
 153 // Use Case:
 154 //  root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1(<pointer_to_phos_geometry_onject>)
 155 //  root [1] cl->Digits2Clusters(digitsTree,clusterTree)
 156 //               //finds RecPoints in the current event
 157 //  root [2] cl->SetDigitsBranch("digits2")
 158 //               //sets another title for Digitis (input) branch
 159 //  root [3] cl->SetRecPointsBranch("recp2")
 160 //               //sets another title four output branches
 161 //  root [4] cl->SetEmcLocalMaxCut(0.03)
 162 //               //set clusterization parameters
 163
 164 // --- ROOT system ---
 165
 166 #include "TMath.h"
 167 #include "TMinuit.h"
 168 #include "TTree.h"
 169 #include "TBenchmark.h"
 170 #include "TClonesArray.h"
 171
 172 // --- Standard library ---
 173
 174 // --- AliRoot header files ---
 175 #include "AliLog.h"
 176 #include "AliConfig.h"
 177 #include "AliPHOSGeometry.h"
 178 #include "AliPHOSClusterizerv1.h"
 179 #include "AliPHOSEmcRecPoint.h"
 180 #include "AliPHOSCpvRecPoint.h"
 181 #include "AliPHOSDigit.h"
 182 #include "AliPHOSDigitizer.h"
 183 #include "AliPHOSCalibrationDB.h"
 184 #include "AliCDBManager.h"
 185 #include "AliCDBStorage.h"
 186 #include "AliCDBEntry.h"
 187 #include "AliPHOSRecoParam.h"
 188 #include "AliPHOSReconstructor.h"
 189 #include "AliPHOSCalibData.h"
 190
 191 ClassImp(AliPHOSClusterizerv1)
 192
 193 //____________________________________________________________________________
 194 AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
 195   AliPHOSClusterizer(),
 196   fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
 197   fWrite(0),
 198   fNumberOfEmcClusters(0),    fNumberOfCpvClusters(0),
 199   fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
 200   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
 201   fW0CPV(0),                  fEmcTimeGate(0),          fEcoreRadius(0)
 202 {
 203   // default ctor (to be used mainly by Streamer)
 204
 205   fDefaultInit = kTRUE ;
 206 }
 207
 208 //____________________________________________________________________________
 209 AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
 210   AliPHOSClusterizer(geom),
 211   fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
 212   fWrite(0),
 213   fNumberOfEmcClusters(0),    fNumberOfCpvClusters(0),
 214   fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
 215   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
 216   fW0CPV(0),                  fEmcTimeGate(0),          fEcoreRadius(0)
 217 {
 218   // ctor with the indication of the file where header Tree and digits Tree are stored
 219
 220   Init() ;
 221   fDefaultInit = kFALSE ;
 222 }
 223
 224 //____________________________________________________________________________
 225   AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
 226 {
 227   // dtor
 228
 229 }
 230 //____________________________________________________________________________
 231 void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
 232 {
 233   // Steering method to perform clusterization for one event
 234   // The input is the tree with digits.
 235   // The output is the tree with clusters.
 236
 237   if(strstr(option,"tim"))
 238     gBenchmark->Start("PHOSClusterizer");
 239
 240   if(strstr(option,"print")) {
 241     Print() ;
 242     return ;
 243   }
 244
 245   MakeClusters() ;
 246
 247   AliDebug(2,Form(" ---- Printing clusters (%d)\n",
 248                   fEMCRecPoints->GetEntries()));
 249   if(AliLog::GetGlobalDebugLevel()>1)
 250     fEMCRecPoints->Print();
 251
 252   if(fToUnfold)
 253     MakeUnfolding();
 254
 255   WriteRecPoints();
 256
 257   if(strstr(option,"deb"))
 258     PrintRecPoints(option) ;
 259
 260   if(strstr(option,"tim")){
 261     gBenchmark->Stop("PHOSClusterizer");
 262     AliInfo(Form("took %f seconds for Clusterizing\n",
 263                  gBenchmark->GetCpuTime("PHOSClusterizer")));
 264   }
 265   fEMCRecPoints->Delete();
 266   fCPVRecPoints->Delete();
 267 }
 268
 269 //____________________________________________________________________________
 270 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
 271                                     Int_t nPar, Float_t * fitparameters) const
 272 {
 273   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 274   // The initial values for fitting procedure are set equal to the positions of local maxima.
 275   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 276
 277
 278   gMinuit->mncler();                     // Reset Minuit's list of paramters
 279   gMinuit->SetPrintLevel(-1) ;           // No Printout
 280   gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
 281                                          // To set the address of the minimization function
 282
 283   TList * toMinuit = new TList();
 284   toMinuit->AddAt(emcRP,0) ;
 285   toMinuit->AddAt(fDigitsArr,1) ;
 286   toMinuit->AddAt(fGeom,2) ;
 287
 288   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 289
 290   // filling initial values for fit parameters
 291   AliPHOSDigit * digit ;
 292
 293   Int_t ierflg  = 0;
 294   Int_t index   = 0 ;
 295   Int_t nDigits = (Int_t) nPar / 3 ;
 296
 297   Int_t iDigit ;
 298
 299   for(iDigit = 0; iDigit < nDigits; iDigit++){
 300     digit = maxAt[iDigit];
 301
 302     Int_t relid[4] ;
 303     Float_t x = 0.;
 304     Float_t z = 0.;
 305     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 306     fGeom->RelPosInModule(relid, x, z) ;
 307
 308     Float_t energy = maxAtEnergy[iDigit] ;
 309
 310     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 311     index++ ;
 312     if(ierflg != 0){
 313       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
 314       return kFALSE;
 315     }
 316     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 317     index++ ;
 318     if(ierflg != 0){
 319        Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
 320       return kFALSE;
 321     }
 322     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 323     index++ ;
 324     if(ierflg != 0){
 325       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
 326       return kFALSE;
 327     }
 328   }
 329
 330   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 331                       //  depends on it.
 332   Double_t p1 = 1.0 ;
 333   Double_t p2 = 0.0 ;
 334
 335   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 336   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 337   gMinuit->SetMaxIterations(5);
 338   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 339
 340   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 341
 342   if(ierflg == 4){  // Minimum not found
 343     Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
 344     return kFALSE ;
 345   }
 346   for(index = 0; index < nPar; index++){
 347     Double_t err ;
 348     Double_t val ;
 349     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 350     fitparameters[index] = val ;
 351    }
 352
 353   delete toMinuit ;
 354   return kTRUE;
 355
 356 }
 357
 358
 359 //____________________________________________________________________________
 360 void AliPHOSClusterizerv1::Init()
 361 {
 362   // Make all memory allocations which can not be done in default constructor.
 363   // Attach the Clusterizer task to the list of PHOS tasks
 364
 365   fEmcCrystals = fGeom->GetNModules() *  fGeom->GetNCristalsInModule() ;
 366
 367   if(!gMinuit)
 368     gMinuit = new TMinuit(100);
 369
 370   if (!fgCalibData)
 371     fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
 372   if (fgCalibData->GetCalibDataEmc() == 0)
 373     AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
 374   if (fgCalibData->GetCalibDataCpv() == 0)
 375     AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
 376
 377 }
 378
 379 //____________________________________________________________________________
 380 void AliPHOSClusterizerv1::InitParameters()
 381 {
 382
 383   fNumberOfCpvClusters     = 0 ;
 384   fNumberOfEmcClusters     = 0 ;
 385
 386   const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
 387   if(!recoParam) AliFatal("Reconstruction parameters are not set!");
 388
 389   recoParam->Print();
 390
 391   fEmcClusteringThreshold  = recoParam->GetEMCClusteringThreshold();
 392   fCpvClusteringThreshold  = recoParam->GetCPVClusteringThreshold();
 393
 394   fEmcLocMaxCut            = recoParam->GetEMCLocalMaxCut();
 395   fCpvLocMaxCut            = recoParam->GetCPVLocalMaxCut();
 396
 397   fW0                      = recoParam->GetEMCLogWeight();
 398   fW0CPV                   = recoParam->GetCPVLogWeight();
 399
 400   fEmcTimeGate             = 1.e-6 ;
 401   fEcoreRadius             = recoParam->GetEMCEcoreRadius();
 402
 403   fToUnfold                = recoParam->EMCToUnfold() ;
 404
 405   fWrite                   = kTRUE ;
 406 }
 407
 408 //____________________________________________________________________________
 409 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
 410 {
 411   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 412   //                                       = 1 are neighbour
 413   //                                       = 2 are not neighbour but do not continue searching
 414   //                                       =-1 are not neighbour, continue searching, but do not look before d2 next time
 415   // neighbours are defined as digits having at least a common vertex
 416   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 417   //                                      which is compared to a digit (d2)  not yet in a cluster
 418
 419   Int_t relid1[4] ;
 420   fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
 421
 422   Int_t relid2[4] ;
 423   fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
 424
 425   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
 426     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 427     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 428
 429     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){   //At least common vertex
 430       //    if (( relid1[2]==relid2[2] && coldiff <= 1 )  || ( relid1[3]==relid2[3] &&  rowdiff <= 1 )){ //common side
 431       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
 432       return 1 ;
 433     }
 434     else {
 435       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
 436         return 2; //  Difference in row numbers is too large to look further
 437     }
 438     return 0 ;
 439
 440   }
 441   else {
 442     if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
 443       return -1 ;
 444     if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
 445       return -1 ;
 446
 447     return 2 ;
 448
 449   }
 450
 451   return 0 ;
 452 }
 453 //____________________________________________________________________________
 454 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
 455 {
 456   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
 457
 458   Bool_t rv = kFALSE ;
 459
 460   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 461
 462   if(digit->GetId() <= nEMC )   rv = kTRUE;
 463
 464   return rv ;
 465 }
 466
 467 //____________________________________________________________________________
 468 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 469 {
 470   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
 471
 472   Bool_t rv = kFALSE ;
 473
 474   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 475
 476   if(digit->GetId() > nEMC )   rv = kTRUE;
 477
 478   return rv ;
 479 }
 480
 481 //____________________________________________________________________________
 482 void AliPHOSClusterizerv1::WriteRecPoints()
 483 {
 484
 485   // Creates new branches with given title
 486   // fills and writes into TreeR.
 487
 488   Int_t index ;
 489   //Evaluate position, dispersion and other RecPoint properties..
 490   Int_t nEmc = fEMCRecPoints->GetEntriesFast();
 491   Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
 492   TVector3 fakeVtx(0.,0.,0.) ;
 493   for(index = 0; index < nEmc; index++){
 494     AliPHOSEmcRecPoint * rp =
 495       dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
 496     rp->Purify(emcMinE) ;
 497     if(rp->GetMultiplicity()==0){
 498       fEMCRecPoints->RemoveAt(index) ;
 499       delete rp ;
 500       continue;
 501     }
 502
 503     // No vertex is available now, calculate corrections in PID
 504     rp->EvalAll(fDigitsArr) ;
 505     rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
 506     rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
 507     rp->EvalLocal2TrackingCSTransform();
 508   }
 509   fEMCRecPoints->Compress() ;
 510   fEMCRecPoints->Sort() ;
 511   //  fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
 512   for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
 513     dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
 514   }
 515
 516   //For each rec.point set the distance to the nearest bad crystal (BVP)
 517   SetDistancesToBadChannels();
 518
 519   //Now the same for CPV
 520   for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
 521     AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
 522     rp->EvalAll(fDigitsArr) ;
 523     rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
 524     rp->EvalLocal2TrackingCSTransform();
 525   }
 526   fCPVRecPoints->Sort() ;
 527
 528   for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
 529     dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
 530
 531   fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
 532
 533   if(fWrite){ //We write TreeR
 534     fTreeR->Fill();
 535   }
 536 }
 537
 538 //____________________________________________________________________________
 539 void AliPHOSClusterizerv1::MakeClusters()
 540 {
 541   // Steering method to construct the clusters stored in a list of Reconstructed Points
 542   // A cluster is defined as a list of neighbour digits
 543
 544   fNumberOfCpvClusters     = 0 ;
 545   fNumberOfEmcClusters     = 0 ;
 546
 547   //Mark all digits as unused yet
 548   const Int_t maxNDigits = 3584; // There is no clusters larger than PHOS module ;)
 549   Int_t nDigits=fDigitsArr->GetEntriesFast() ;
 550
 551   for(Int_t i=0; i<nDigits; i++){
 552     fDigitsUsed[i]=0 ;
 553   }
 554   Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
 555                      //e.g. first digit in this module, first CPV digit etc.
 556   AliPHOSDigit * digit ;
 557   TArrayI clusterdigitslist(maxNDigits) ;
 558   AliPHOSRecPoint * clu = 0 ;
 559   for(Int_t i=0; i<nDigits; i++){
 560     if(fDigitsUsed[i])
 561       continue ;
 562
 563     digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
 564
 565     clu=0 ;
 566
 567     Int_t index ;
 568
 569     //is this digit so energetic that start cluster?
 570     if (( IsInEmc(digit) &&  Calibrate(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
 571         ( IsInCpv(digit) &&  Calibrate(digit->GetEnergy(),digit->GetId()) > fCpvClusteringThreshold ) ) {
 572       Int_t iDigitInCluster = 0 ;
 573       if  ( IsInEmc(digit) ) {
 574         // start a new EMC RecPoint
 575         if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 576           fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
 577
 578         fEMCRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
 579         clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 580         fNumberOfEmcClusters++ ;
 581         clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId())) ;
 582         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
 583         iDigitInCluster++ ;
 584         fDigitsUsed[i]=kTRUE ;
 585       } else {
 586         // start a new CPV cluster
 587         if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 588           fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
 589
 590         fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
 591         clu =  static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 592         fNumberOfCpvClusters++ ;
 593         clu->AddDigit(*digit,  Calibrate(digit->GetEnergy(),digit->GetId())) ;
 594         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;
 595         iDigitInCluster++ ;
 596         fDigitsUsed[i]=kTRUE ;
 597       } // else
 598
 599       //Now scan remaining digits in list to find neigbours of our seed
 600
 601       AliPHOSDigit * digitN ;
 602       index = 0 ;
 603       while (index < iDigitInCluster){ // scan over digits already in cluster
 604         digit =  static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) )  ;
 605         index++ ;
 606         for(Int_t j=iFirst; j<nDigits; j++){
 607           if (iDigitInCluster >= maxNDigits) {
 608             AliError(Form("The number of digits in cluster is more than %d, skip the rest of event",
 609                           maxNDigits));
 610             return;
 611           }
 612           if(fDigitsUsed[j])
 613             continue ;        //look through remaining digits
 614           digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
 615           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 616           switch (ineb ) {
 617           case -1:   //too early (e.g. previous module), do not look before j at subsequent passes
 618             iFirst=j ;
 619             break ;
 620           case 0 :   // not a neighbour
 621             break ;
 622           case 1 :   // are neighbours
 623             clu->AddDigit(*digitN, Calibrate(digitN->GetEnergy(),digit->GetId())) ;
 624             clusterdigitslist[iDigitInCluster] = j ;
 625             iDigitInCluster++ ;
 626             fDigitsUsed[j]=kTRUE ;
 627             break ;
 628           case 2 :   // too far from each other
 629             goto endOfLoop;
 630           } // switch
 631
 632         }
 633
 634         endOfLoop: ; //scanned all possible neighbours for this digit
 635
 636       } // loop over cluster
 637     } // energy theshold
 638   }
 639
 640 }
 641
 642 //____________________________________________________________________________
 643 void AliPHOSClusterizerv1::MakeUnfolding()
 644 {
 645   // Unfolds clusters using the shape of an ElectroMagnetic shower
 646   // Performs unfolding of all EMC/CPV clusters
 647
 648   // Unfold first EMC clusters
 649   if(fNumberOfEmcClusters > 0){
 650
 651     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 652
 653     Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
 654     Int_t index ;
 655     for(index = 0 ; index < numberofNotUnfolded ; index++){
 656
 657       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
 658       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
 659         break ;
 660
 661       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 662       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 663       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 664       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
 665
 666       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 667         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 668
 669         fEMCRecPoints->Remove(emcRecPoint);
 670         fEMCRecPoints->Compress() ;
 671         index-- ;
 672         fNumberOfEmcClusters -- ;
 673         numberofNotUnfolded-- ;
 674       }
 675       else{
 676         emcRecPoint->SetNExMax(1) ; //Only one local maximum
 677       }
 678
 679       delete[] maxAt ;
 680       delete[] maxAtEnergy ;
 681     }
 682   }
 683   // Unfolding of EMC clusters finished
 684
 685
 686   // Unfold now CPV clusters
 687   if(fNumberOfCpvClusters > 0){
 688
 689     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 690
 691     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
 692     Int_t index ;
 693     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
 694
 695       AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
 696
 697       if(recPoint->GetPHOSMod()> nModulesToUnfold)
 698         break ;
 699
 700       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
 701
 702       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 703       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 704       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 705       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
 706
 707       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 708         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 709         fCPVRecPoints->Remove(emcRecPoint);
 710         fCPVRecPoints->Compress() ;
 711         index-- ;
 712         numberofCpvNotUnfolded-- ;
 713         fNumberOfCpvClusters-- ;
 714       }
 715
 716       delete[] maxAt ;
 717       delete[] maxAtEnergy ;
 718     }
 719   }
 720   //Unfolding of Cpv clusters finished
 721
 722 }
 723
 724 //____________________________________________________________________________
 725 Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
 726 {
 727   // Shape of the shower (see PHOS TDR)
 728   // If you change this function, change also the gradient evaluation in ChiSquare()
 729
 730   //for the moment we neglect dependence on the incident angle.
 731
 732   Double_t r2    = x*x + z*z ;
 733   Double_t r4    = r2*r2 ;
 734   Double_t r295  = TMath::Power(r2, 2.95/2.) ;
 735   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 736   return shape ;
 737 }
 738
 739 //____________________________________________________________________________
 740 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 741                                           Int_t nMax,
 742                                           AliPHOSDigit ** maxAt,
 743                                           Float_t * maxAtEnergy)
 744 {
 745   // Performs the unfolding of a cluster with nMax overlapping showers
 746
 747   Int_t nPar = 3 * nMax ;
 748   Float_t * fitparameters = new Float_t[nPar] ;
 749
 750   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
 751
 752   if( !rv ) {
 753     // Fit failed, return and remove cluster
 754     iniEmc->SetNExMax(-1) ;
 755     delete[] fitparameters ;
 756     return ;
 757   }
 758
 759   // create ufolded rec points and fill them with new energy lists
 760   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
 761   // and later correct this number in acordance with actual energy deposition
 762
 763   Int_t nDigits = iniEmc->GetMultiplicity() ;
 764   Float_t * efit = new Float_t[nDigits] ;
 765   Float_t xDigit=0.,zDigit=0. ;
 766   Float_t xpar=0.,zpar=0.,epar=0.  ;
 767   Int_t relid[4] ;
 768   AliPHOSDigit * digit = 0 ;
 769   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 770
 771   TVector3 vIncid ;
 772
 773   Int_t iparam ;
 774   Int_t iDigit ;
 775   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 776     digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
 777     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 778     fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 779     efit[iDigit] = 0;
 780
 781     iparam = 0 ;
 782     while(iparam < nPar ){
 783       xpar = fitparameters[iparam] ;
 784       zpar = fitparameters[iparam+1] ;
 785       epar = fitparameters[iparam+2] ;
 786       iparam += 3 ;
 787 //      fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
 788 //      efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
 789       efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 790     }
 791   }
 792
 793   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
 794   // so that energy deposited in each cell is distributed betwin new clusters proportionally
 795   // to its contribution to efit
 796
 797   Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
 798   Float_t ratio ;
 799
 800   iparam = 0 ;
 801   while(iparam < nPar ){
 802     xpar = fitparameters[iparam] ;
 803     zpar = fitparameters[iparam+1] ;
 804     epar = fitparameters[iparam+2] ;
 805     iparam += 3 ;
 806 //    fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
 807
 808     AliPHOSEmcRecPoint * emcRP = 0 ;
 809
 810     if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
 811
 812       if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 813         fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
 814
 815       (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
 816       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 817       fNumberOfEmcClusters++ ;
 818       emcRP->SetNExMax((Int_t)nPar/3) ;
 819     }
 820     else{//create new entries in fCPVRecPoints
 821       if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 822         fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
 823
 824       (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
 825       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 826       fNumberOfCpvClusters++ ;
 827     }
 828
 829     Float_t eDigit ;
 830     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 831       digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
 832       fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 833       fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 834 //      ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
 835       ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
 836       eDigit = emcEnergies[iDigit] * ratio ;
 837       emcRP->AddDigit( *digit, eDigit ) ;
 838     }
 839   }
 840
 841   delete[] fitparameters ;
 842   delete[] efit ;
 843
 844 }
 845
 846 //_____________________________________________________________________________
 847 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 848 {
 849   // Calculates the Chi square for the cluster unfolding minimization
 850   // Number of parameters, Gradient, Chi squared, parameters, what to do
 851
 852   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
 853
 854   AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
 855   TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
 856   // A bit buggy way to get an access to the geometry
 857   // To be revised!
 858   AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
 859
 860 //  TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ;  //Vertex position
 861
 862   //  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
 863
 864   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 865
 866   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
 867
 868   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
 869
 870 //  TVector3 vInc ;
 871   fret = 0. ;
 872   Int_t iparam ;
 873
 874   if(iflag == 2)
 875     for(iparam = 0 ; iparam < nPar ; iparam++)
 876       Grad[iparam] = 0 ; // Will evaluate gradient
 877
 878   Double_t efit ;
 879
 880   AliPHOSDigit * digit ;
 881   Int_t iDigit ;
 882
 883   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 884
 885     digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
 886
 887     Int_t relid[4] ;
 888     Float_t xDigit ;
 889     Float_t zDigit ;
 890
 891     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 892
 893     geom->RelPosInModule(relid, xDigit, zDigit) ;
 894
 895      if(iflag == 2){  // calculate gradient
 896        Int_t iParam = 0 ;
 897        efit = 0 ;
 898        while(iParam < nPar ){
 899          Double_t dx = (xDigit - x[iParam]) ;
 900          iParam++ ;
 901          Double_t dz = (zDigit - x[iParam]) ;
 902          iParam++ ;
 903 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
 904 //         efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
 905          efit += x[iParam] * ShowerShape(dx,dz) ;
 906          iParam++ ;
 907        }
 908        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
 909        iParam = 0 ;
 910        while(iParam < nPar ){
 911          Double_t xpar = x[iParam] ;
 912          Double_t zpar = x[iParam+1] ;
 913          Double_t epar = x[iParam+2] ;
 914 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
 915          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
 916 //         Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
 917          Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 918 //DP: No incident angle dependence in gradient yet!!!!!!
 919          Double_t r4 = dr*dr*dr*dr ;
 920          Double_t r295 = TMath::Power(dr,2.95) ;
 921          Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
 922                                          0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
 923
 924          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
 925          iParam++ ;
 926          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
 927          iParam++ ;
 928          Grad[iParam] += shape ;                                  // Derivative over energy
 929          iParam++ ;
 930        }
 931      }
 932      efit = 0;
 933      iparam = 0 ;
 934
 935      while(iparam < nPar ){
 936        Double_t xpar = x[iparam] ;
 937        Double_t zpar = x[iparam+1] ;
 938        Double_t epar = x[iparam+2] ;
 939        iparam += 3 ;
 940 //       fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
 941 //       efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
 942        efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 943      }
 944
 945      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
 946      // Here we assume, that sigma = sqrt(E)
 947   }
 948
 949 }
 950
 951 //____________________________________________________________________________
 952 void AliPHOSClusterizerv1::Print(const Option_t *)const
 953 {
 954   // Print clusterizer parameters
 955
 956   TString message ;
 957   TString taskName(GetName()) ;
 958   taskName.ReplaceAll(Version(), "") ;
 959
 960   if( strcmp(GetName(), "") !=0 ) {
 961     // Print parameters
 962     message  = "\n--------------- %s %s -----------\n" ;
 963     message += "Clusterizing digits from the file: %s\n" ;
 964     message += "                           Branch: %s\n" ;
 965     message += "                       EMC Clustering threshold = %f\n" ;
 966     message += "                       EMC Local Maximum cut    = %f\n" ;
 967     message += "                       EMC Logarothmic weight   = %f\n" ;
 968     message += "                       CPV Clustering threshold = %f\n" ;
 969     message += "                       CPV Local Maximum cut    = %f\n" ;
 970     message += "                       CPV Logarothmic weight   = %f\n" ;
 971     if(fToUnfold)
 972       message += " Unfolding on\n" ;
 973     else
 974       message += " Unfolding off\n" ;
 975
 976     message += "------------------------------------------------------------------" ;
 977   }
 978   else
 979     message = " AliPHOSClusterizerv1 not initialized " ;
 980
 981   AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
 982        taskName.Data(),
 983        GetTitle(),
 984        taskName.Data(),
 985        GetName(),
 986        fEmcClusteringThreshold,
 987        fEmcLocMaxCut,
 988        fW0,
 989        fCpvClusteringThreshold,
 990        fCpvLocMaxCut,
 991        fW0CPV )) ;
 992 }
 993 //____________________________________________________________________________
 994 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
 995 {
 996   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
 997
 998   AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
 999                fEMCRecPoints->GetEntriesFast(),
1000                fCPVRecPoints->GetEntriesFast() ))  ;
1001
1002   if(strstr(option,"all")) {
1003     printf("\n EMC clusters \n") ;
1004     printf("Index    Ene(MeV) Multi Module     X    Y   Z    Lambdas_1  Lambda_2  # of prim  Primaries list\n") ;
1005     Int_t index ;
1006     for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1007       AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1008       TVector3  locpos;
1009       rp->GetLocalPosition(locpos);
1010       Float_t lambda[2];
1011       rp->GetElipsAxis(lambda);
1012       Int_t * primaries;
1013       Int_t nprimaries;
1014       primaries = rp->GetPrimaries(nprimaries);
1015       printf("\n%6d  %8.2f  %3d     %2d     %4.1f   %4.1f   %4.1f   %4f  %4f    %2d     : ",
1016               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1017               locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1018
1019       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1020         printf("%d ", primaries[iprimary] ) ;
1021       }
1022       printf("\n") ;
1023     }
1024
1025     //Now plot CPV recPoints
1026     printf("\n CPV clusters \n") ;
1027     printf("Index    Ene(MeV) Module     X     Y    Z  \n") ;
1028     for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1029       AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1030
1031       TVector3  locpos;
1032       rp->GetLocalPosition(locpos);
1033
1034       printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f \n",
1035              rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1036              locpos.X(), locpos.Y(), locpos.Z()) ;
1037     }
1038   }
1039 }
1040
1041
1042 //____________________________________________________________________________
1043 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1044 {
1045   //For each EMC rec. point set the distance to the nearest bad crystal.
1046   //Author: Boris Polichtchouk
1047
1048   if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1049   AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
1050
1051   Int_t badIds[8000];
1052   fgCalibData->EmcBadChannelIds(badIds);
1053
1054   AliPHOSEmcRecPoint* rp;
1055
1056   TMatrixF gmat;
1057   TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1058   TVector3 gposBadChannel; // global position of bad crystal
1059   TVector3 dR;
1060
1061   Float_t dist,minDist;
1062   Int_t relid[4]={0,0,0,0} ;
1063   TVector3 lpos ;
1064   for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1065     rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1066     //evaluate distance to border
1067     relid[0]=rp->GetPHOSMod() ;
1068     relid[2]=1 ;
1069     relid[3]=1 ;
1070     Float_t xcorner,zcorner;
1071     fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1072     rp->GetLocalPosition(lpos) ;
1073     minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
1074     for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1075       fGeom->AbsToRelNumbering(badIds[iBad],relid)  ;
1076       if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1077         continue ;                   //bad channels can be in the module which does not exist in simulations.
1078       rp->GetGlobalPosition(gposRecPoint,gmat);
1079       fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1080       AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1081                       gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1082                       gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1083       dR = gposBadChannel-gposRecPoint;
1084       dist = dR.Mag();
1085       if(dist<minDist) minDist = dist;
1086     }
1087
1088     rp->SetDistanceToBadCrystal(minDist);
1089   }
1090
1091 }
1092 //==================================================================================
1093 Float_t AliPHOSClusterizerv1::Calibrate(Float_t amp, Int_t absId){
1094   // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB
1095
1096   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1097
1098   //Determine rel.position of the cell absolute ID
1099   Int_t relId[4];
1100   geom->AbsToRelNumbering(absId,relId);
1101   Int_t module=relId[0];
1102   Int_t row   =relId[2];
1103   Int_t column=relId[3];
1104   if(relId[1]){ //CPV
1105     Float_t calibration = fgCalibData->GetADCchannelCpv(module,column,row);
1106     return amp*calibration ;
1107   }
1108   else{ //EMC
1109     Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row);
1110     return amp*calibration ;
1111   }
1112 }