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